Does a cooperative environment affect isolation?

Does allocating cooperative individuals to high-degree nodes prevent isolation?

round = NULL
gameID = NULL 
nodes = NULL
meanDegree = NULL
meanDegreeCoop = NULL #mean degree among those who cooperated
meanDegreeDefect = NULL #mean degree among those who defected
meanDist = NULL #mean shortest path length divided by the maximum possible path length (Nvertices - 1)
maxDist = NULL #max shortest path length divided by the maximum possible path length (Nvertices - 1)
coefGlobal = NULL #global clustering coefficient
percentCoop = NULL #proportion of those who cooperated
percentCoopPopular = NULL #proportion of those who cooperated, among those with degree >=10
coopEnv = NULL 
assort = NULL

for(i in 1:length(g.isolation)){
    round[i] =  mean(V(g.isolation[[i]])$round)
    gameID[i] = V(g.isolation[[i]])$gameID[1]
    nodes[i] = vcount(g.isolation[[i]])
    meanDegree[i] = mean(igraph::degree(g.isolation[[i]]),na.rm=TRUE)
    meanDegreeCoop[i] = mean(igraph::degree(g.isolation[[i]])[V(g.isolation[[i]])$behavior=="C"],na.rm=TRUE)
    meanDegreeDefect[i] = mean(igraph::degree(g.isolation[[i]])[V(g.isolation[[i]])$behavior=="D"],na.rm=TRUE)
    meanDist[i] = mean(distances(g.isolation[[i]])[distances(g.isolation[[i]])!=0 & distances(g.isolation[[i]])!=Inf],na.rm=TRUE) / length(V(g.isolation[[i]])-1)
    maxDist[i] = max(distances(g.isolation[[i]])[distances(g.isolation[[i]])!=0 & distances(g.isolation[[i]])!=Inf],na.rm=TRUE) / length(V(g.isolation[[i]])-1)
    coefGlobal[i] = transitivity(g.isolation[[i]], type="global")
    percentCoop[i] = prop.table(table(V(g.isolation[[i]])$behavior))["C"]
    percentCoopPopular[i] = prop.table(table(V(g.isolation[[i]])$behavior[igraph::degree(g.isolation[[i]])>=6]))["C"]
    coopEnv[i] = sum(ifelse(V(g.isolation[[i]])$behavior=="C",1,0)*igraph::degree(g.isolation[[i]])/sum(igraph::degree(g.isolation[[i]])),na.rm=TRUE)
    assort[i] = assortativity(g.isolation[[i]], V(g.isolation[[i]])$behavior == "C")
}

smallWorld = data.frame(
  round = round,
  gameID = gameID,
  nodes = nodes,
  meanDegree = meanDegree,
  meanDegreeCoop = meanDegreeCoop,
  meanDegreeDefect = meanDegreeDefect,
  diffDegreeCD = meanDegreeCoop - meanDegreeDefect,
  meanDist = meanDist,
  maxDist = maxDist,
  coefGlobal = coefGlobal,
  percentCoop = percentCoop,
  percentCoopPopular = percentCoopPopular,
  coopEnv = coopEnv,
  assort = assort
)


round = NULL
gameID = NULL
visibleWealth = NULL
n=1
for(i in unique(cdata$round)){
  for(m in unique(cdata$gameID)){
    round[n] = as.numeric(i)-1
    gameID[n] = m
    visibleWealth[n] = ifelse(cdata[cdata$round==i & cdata$gameID==m,]$showScore[1]=="true",1,0)
    n=n+1
  }
}

smallWorld.2 = data.frame(
  round = round,
  gameID = gameID,
  visibleWealth = visibleWealth
)

smallWorld = merge(smallWorld, smallWorld.2, by=c("gameID","round"), all.x=TRUE)

smallWorld.initial = subset(smallWorld, smallWorld$round==0)

#percent of people that end up being isolated for each gameID
percentIsolated = NULL
percentIsolated = aggregate(isolated ~ gameID, data=sample.wide, FUN=max, na.rm=TRUE)
smallWorld.initial = merge(smallWorld.initial[c("gameID","meanDegree","meanDegreeCoop","meanDegreeDefect","diffDegreeCD","meanDist","coefGlobal","visibleWealth","percentCoop","percentCoopPopular","coopEnv","assort")],percentIsolated, by="gameID", all.x=TRUE)

percentIsolated = NULL
percentIsolated = aggregate(isolated ~ gameID, data=sample.wide, FUN=mean, na.rm=TRUE)
percentIsolated = percentIsolated %>% rename(
  percentIsolated = isolated
)
smallWorld.initial = merge(smallWorld.initial,percentIsolated, by="gameID", all.x=TRUE)




reg = glm(isolated ~ coopEnv, 
      data = smallWorld.initial,
      family = binomial(link = "logit"))
summary(reg)

## 
## Call:
## glm(formula = isolated ~ coopEnv, family = binomial(link = "logit"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3756   0.2641   0.4839   0.6926   1.3696  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    6.047      1.734   3.487 0.000489 ***
## coopEnv       -6.813      2.351  -2.898 0.003750 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 82.760  on 79  degrees of freedom
## Residual deviance: 72.371  on 78  degrees of freedom
## AIC: 76.371
## 
## Number of Fisher Scoring iterations: 5

reg = glm(isolated ~ coopEnv*percentCoop, 
      data = smallWorld.initial,
      family = binomial(link = "logit"))
summary(reg)

## 
## Call:
## glm(formula = isolated ~ coopEnv * percentCoop, family = binomial(link = "logit"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2132   0.4220   0.4658   0.5862   1.7551  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)           -3.096      6.866  -0.451    0.652
## coopEnv               10.481     12.639   0.829    0.407
## percentCoop           10.681     12.929   0.826    0.409
## coopEnv:percentCoop  -20.525     14.822  -1.385    0.166
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 82.760  on 79  degrees of freedom
## Residual deviance: 70.472  on 76  degrees of freedom
## AIC: 78.472
## 
## Number of Fisher Scoring iterations: 4

reg = glm(isolated ~ assort, 
      data = smallWorld.initial,
      family = binomial(link = "logit"))
summary(reg)

## 
## Call:
## glm(formula = isolated ~ assort, family = binomial(link = "logit"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9579   0.5132   0.6439   0.7269   1.0178  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1349     0.3207   3.539 0.000402 ***
## assort       -3.7358     3.0781  -1.214 0.224862    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 67.731  on 64  degrees of freedom
## Residual deviance: 66.177  on 63  degrees of freedom
##   (15 observations deleted due to missingness)
## AIC: 70.177
## 
## Number of Fisher Scoring iterations: 4

reg = glm(isolated ~ assort*percentCoop, 
      data = smallWorld.initial,
      family = binomial(link = "logit"))
summary(reg)

## 
## Call:
## glm(formula = isolated ~ assort * percentCoop, family = binomial(link = "logit"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5706   0.1873   0.4272   0.6978   1.3591  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)   
## (Intercept)           6.695      2.392   2.799  0.00512 **
## assort              -23.388     22.419  -1.043  0.29684   
## percentCoop          -7.680      3.186  -2.411  0.01593 * 
## assort:percentCoop   25.842     29.451   0.877  0.38024   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 67.731  on 64  degrees of freedom
## Residual deviance: 56.386  on 61  degrees of freedom
##   (15 observations deleted due to missingness)
## AIC: 64.386
## 
## Number of Fisher Scoring iterations: 5

reg = glm(percentIsolated ~ coopEnv, 
      data = smallWorld.initial, 
      family = gaussian(link = "identity"))
summary(reg)

## 
## Call:
## glm(formula = percentIsolated ~ coopEnv, family = gaussian(link = "identity"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.12622  -0.05564  -0.02541   0.03834   0.25903  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.20852    0.04303   4.846 6.29e-06 ***
## coopEnv     -0.17058    0.06391  -2.669  0.00925 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.006957149)
## 
##     Null deviance: 0.59222  on 79  degrees of freedom
## Residual deviance: 0.54266  on 78  degrees of freedom
## AIC: -166.43
## 
## Number of Fisher Scoring iterations: 2

reg = glm(percentIsolated ~ coopEnv*percentCoop, 
      data = smallWorld.initial, 
      family = gaussian(link = "identity"))
summary(reg)

## 
## Call:
## glm(formula = percentIsolated ~ coopEnv * percentCoop, family = gaussian(link = "identity"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.15277  -0.04970  -0.01807   0.03678   0.25172  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)           0.1084     0.1925   0.563    0.575
## coopEnv               0.4950     0.3742   1.323    0.190
## percentCoop          -0.2223     0.3692  -0.602    0.549
## coopEnv:percentCoop  -0.3977     0.4417  -0.900    0.371
## 
## (Dispersion parameter for gaussian family taken to be 0.006743275)
## 
##     Null deviance: 0.59222  on 79  degrees of freedom
## Residual deviance: 0.51249  on 76  degrees of freedom
## AIC: -167.01
## 
## Number of Fisher Scoring iterations: 2

reg = glm(percentIsolated ~ assort, 
      data = smallWorld.initial, 
      family = gaussian(link = "identity"))
summary(reg)

## 
## Call:
## glm(formula = percentIsolated ~ assort, family = gaussian(link = "identity"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.13800  -0.04889  -0.01599   0.03676   0.29844  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.08136    0.01164   6.988  2.1e-09 ***
## assort      -0.33980    0.09891  -3.435  0.00105 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.006992089)
## 
##     Null deviance: 0.52302  on 64  degrees of freedom
## Residual deviance: 0.44050  on 63  degrees of freedom
##   (15 observations deleted due to missingness)
## AIC: -134.16
## 
## Number of Fisher Scoring iterations: 2

reg = glm(percentIsolated ~ assort*percentCoop, 
      data = smallWorld.initial, 
      family = gaussian(link = "identity"))
summary(reg)

## 
## Call:
## glm(formula = percentIsolated ~ assort * percentCoop, family = gaussian(link = "identity"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.13441  -0.04183  -0.01733   0.02556   0.27510  
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.22289    0.05970   3.733 0.000418 ***
## assort             -0.70774    0.42639  -1.660 0.102080    
## percentCoop        -0.20362    0.08504  -2.394 0.019731 *  
## assort:percentCoop  0.59881    0.63264   0.947 0.347615    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.006108948)
## 
##     Null deviance: 0.52302  on 64  degrees of freedom
## Residual deviance: 0.37265  on 61  degrees of freedom
##   (15 observations deleted due to missingness)
## AIC: -141.04
## 
## Number of Fisher Scoring iterations: 2

# #MSM
# # estimation of denominator of ip weights
# #cooperatorの数以外は全てRandom（Erdos-Renyi random graphにより決まるため）なので、coopEnvはpercentCoopで調整すれば十分か
# den.fit.obj <- lm(coopEnv ~ percentCoop + I(percentCoop^2), data = smallWorld.initial)
# p.den <- predict(den.fit.obj, type = "response")
# dens.den <- dnorm(smallWorld.initial$coopEnv, p.den, summary(den.fit.obj)$sigma)
# 
# # estimation of numerator of ip weights
# num.fit.obj <- lm(coopEnv ~ 1, data = smallWorld.initial)
# p.num <- predict(num.fit.obj, type = "response")
# dens.num <- dnorm(smallWorld.initial$coopEnv, p.num, summary(num.fit.obj)$sigma)
# 
# smallWorld.initial$sw.a = dens.num/dens.den
# summary(smallWorld.initial$sw.a)
# 
# msm.sw.cont <-geepack::geeglm(isolated ~ coopEnv, 
#                  data=smallWorld.initial, weights=sw.a, family = binomial(link="logit"), id=gameID, corstr="independence")
# summary(msm.sw.cont)
# 
# msm.sw.cont <-geepack::geeglm(percentIsolated ~ coopEnv, 
#                  data=smallWorld.initial, weights=sw.a, family = gaussian(link="identity"), id=gameID, corstr="independence")
# summary(msm.sw.cont)


ggplot(data = smallWorld.initial, aes(x = coopEnv, y = percentIsolated)) + 
  geom_point() +
  ggtitle("Percent isolated") +
  scale_x_continuous("Cooporative environment") +  
  scale_y_continuous("Percent of individuals isolated") + 
  scale_fill_continuous(type = "viridis")

ggplot(data = smallWorld.initial, aes(x = assort, y = percentIsolated)) + 
  geom_point() +
  ggtitle("Percent isolated") +
  scale_x_continuous("Assortativity coefficient of initial C/D") +  
  scale_y_continuous("Percent of individuals isolated") + 
  scale_fill_continuous(type = "viridis")

## Warning: Removed 15 rows containing missing values (geom_point).

Simulation 1

All graphs start as Erdos-Renyi random newtorks

Change where to place the cooperators in the network: 1. cooperators have low degrees

#Define degrees of isolation
isolationDegree = 2

#number of iterations per arm
iterations = 500


df.netIntLowDegree = data.frame(
  coopFrac = NULL,
  avgCoop = NULL, 
  avgCoopFinal = NULL, 
  percentIsolation = NULL,
  isolation = NULL,
  percentIsolatedD = NULL,
  nCommunities = NULL,
  communitySize = NULL,
  assortativityInitial = NULL,
  assortativityFinal = NULL,
  conversionRate = NULL,
  conversionToD = NULL, 
  conversionToC = NULL, 
  degreeC = NULL, 
  degreeD = NULL,
  meanConversionToD = NULL, 
  meanConversionToC = NULL, 
  degreeLost = NULL,
  degreeLostC = NULL,
  degreeLostD = NULL
)


for(frac in c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)){
  #nodes in the top fractionCoop degrees will automatically be a cooperator
  fractionCoop = frac
  
    coopFrac = NULL
    avgCoop = NULL
    avgCoopFinal = NULL
    percentIsolation = NULL
    isolation = NULL
    percentIsolatedD = NULL
    nCommunities = NULL
    communitySize = NULL
    assortativityInitial = NULL
    assortativityFinal = NULL
    conversionRate = NULL
    conversionToD = NULL 
    conversionToC = NULL
    degreeC = NULL
    degreeD = NULL
    meanConversionToD = NULL 
    meanConversionToC = NULL
    degreeLost = NULL
    degreeLostC = NULL
    degreeLostD = NULL
    for(m in c(1:iterations)){
      # Section 1. NOTES, packages, and Parameters
      #Importing library
      library(igraph) # for network graphing
      library(reldist) # for gini calculatio
      library(boot) # for inv.logit calculation
      #Two prefixed functions
      #rank
      rank1 = function(x) {rank(x,na.last=NA,ties.method="average")[1]} #a smaller value has a smaller rank.
      #gini mean difference (a.k.a. mean difference: please refer to https://stat.ethz.ch/pipermail/r-help/2003-April/032782.html)
      gmd = function(x) {
       x1 = na.omit(x)
       n = length(x1)
      tmp = 0
       for (i in 1:n) {
       for (j in 1:n) {
       tmp <- tmp + abs(x1[i]-x1[j])
       }
       }
      answer = tmp/(n*n)
       return(answer)
      }
      # List of manipulating parameters of experiments
      #L : number of round
      #V : Visible or not
      #A : Income of a rich-group subject
      #B : Income of a poor-group subject
      #R : Probability to be assigned to a rich group
      #I : Number of the same-parameter trial
      #Example
      L = 10
      V = 0
      A = 700
      B = 300
      R = 0.5
      I = 0
      # List of fixed parameters of experiments (assumptions)
      #Rewiring rate = 0.3
      #GINI coefficient (can be known by A or B)
      GINI = 0*as.numeric(A==500) + 0.2*as.numeric(A %in% c(700,850)) + 0.4*as.numeric(A ==1150)
      #Collecting data frame (final output data frame)
      result = data.frame(round=0:L,n_par=NA,n_A=NA,avg_coop=NA,avg_degree=NA,avg_wealth=NA,gini=NA,gmd=NA,avg_coop_A=NA,avg_degree_A=NA,avg_wealth_A=NA,gini_A=NA,gmd_A=NA,avg_coop_B=NA,avg_degree_B=NA,avg_wealth_B=NA,gini_B=NA,gmd_B=NA,isolation=NA,percentIsolation=NA,meanConversionToD=NA,meanConversionToC=NA,degreeLost=NA,degreeLostC=NA,degreeLostD=NA)
      #_A is for a richer group and _B is for a poorer group
      
      
      #####################################################
      # Section 1.5: Practice rounds 1 to 2, to determine C/D in round 1
      
      N = 17 # median of the number of participants over rounds.
      node_rp0 = data.frame(ego_id=1:N, round=0)
      node_import = node_rp0
      
      for (k in 1:2){
        node_rX = node_import #Importing data
        node_rX$round = node_rX$round + 1
        node_rX[is.na(node_rX$prev_degree)==1,"prev_degree"] = 0
        node_rX[is.na(node_rX$prev_local_rate_coop)==1,"prev_local_rate_coop"] = 0
        #Only this calculation needs to change from Round 1
        if (k==1) {
          node_rX$prob_coop = inv.logit(1.099471) 
        } else {
          node_rX$prob_coop = inv.logit((-0.02339288) + (1.46068980)*as.numeric(node_rX$prev_coop==1))
        } 
        
        node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
        
        node_rX$prev_coop = node_rX$coop
          
        assign(paste("coop_rp",k, sep=""),node_rX$coop)
        
        #For the loop
        node_import = node_rX
      }
      
      #cooperation rate in the practice rounds
      coop_rp = apply(cbind(coop_rp1,coop_rp2),1,mean)

      #####################################################
      # Section 2: Round 0 (Agents and environments)
      #Node data generation
      N = 17 # median of the number of participants over rounds.
      node_r0 = data.frame(ego_id=1:N, round=0)
      node_r0$group = sample(c("rich","poor"),N,replace=TRUE,prob=c(R,1-R)) #R is defined as the probability to be assigned to the rich group
      node_r0$initial_wealth = ifelse(node_r0$group=="rich",A,B)
      #Link data generation
      ego_list = NULL
      for (i in 1:N) { ego_list = c(ego_list,rep(i,N)) }
      link_r0 = data.frame(ego_id=ego_list,alt_id=rep(1:N,N))
      link_r0 = link_r0[(link_r0$ego_id < link_r0$alt_id),] #The link was bidirectional, and thus the half and self are omitted.
      link_r0$connected = sample(0:1,dim(link_r0)[1],replace=TRUE,prob=c(0.7,0.3)) #Initial rewiring rate is fixed, 0.3
      link_r0c_ego = link_r0[link_r0$connected==1,]
      link_r0c_alt = link_r0[link_r0$connected==1,]
      colnames(link_r0c_alt) = c("alt_id","ego_id","connected")
      link_r0c = rbind(link_r0c_ego,link_r0c_alt) #this is bidirectional (double counted) for connected ties.
      link_r0c = link_r0c[order(link_r0c$ego_id),]
      link_r0c$alternumber = NA #putting the number for each alter in the same ego
      link_r0c[1,]$alternumber = 1
      for (i in 1:(dim(link_r0c)[1]-1))
        {if (link_r0c[i,]$ego_id == link_r0c[i+1,]$ego_id)
          {link_r0c[i+1,]$alternumber = link_r0c[i,]$alternumber + 1}
        else
          {link_r0c[i+1,]$alternumber = 1}
        #print(i)
        }
      link_r0c2 = reshape(link_r0c, direction = "wide", idvar=c("ego_id","connected"), timevar="alternumber")
      link_r0c2$initial_degree = apply(link_r0c2[,colnames(link_r0c2)[substr(colnames(link_r0c2),1,6) == "alt_id"]],1,function(x){length(na.omit(x))}) #Degree of each ego
      link_r0c2[is.na(link_r0c2$initial_degree)==1,"initial_degree"] = 0
      #Reflect the degree and initial local gini coefficient into the node data
      node_r0 = merge(x=node_r0,y=link_r0c2,all.x=TRUE,all.y=FALSE,by="ego_id")
      node_r0$initial_avg_env_wealth = NA
      node_r0$initial_local_gini = NA #local gini coefficient of the ego and connecting alters
      node_r0$initial_rel_rank = NA #local rank of ego among the ego and connecting alters (divided by the number of the go and connecting alters)
      for (i in 1:(dim(node_r0)[1])){
        node_r0[i,]$initial_avg_env_wealth = mean(na.omit(node_r0[node_r0$ego_id %in%
        node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6) %in% c("ego_id","alt_id")]],"initial_wealth"]))
        node_r0[i,]$initial_local_gini = gini(na.omit(node_r0[node_r0$ego_id %in% node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6)
        %in% c("ego_id","alt_id")]],"initial_wealth"]))
        node_r0[i,]$initial_rel_rank = rank1(na.omit(node_r0[node_r0$ego_id %in% node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6)
        %in% c("ego_id","alt_id")]],"initial_wealth"]))/length(na.omit(node_r0[node_r0$ego_id %in%
        node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6) %in% c("ego_id","alt_id")]],"initial_wealth"]))
        }
      #Finalization of round 0 and Visualization
      #plot(graph.data.frame(link_r0[link_r0$connected==1,],directed=F)) #plot.igraph
      node_r0$everIsolated = 0
      node_r0$maxDegreeLost = NA
      result[result$round==0,2:25] = c(length(node_r0$ego_id),length(node_r0[node_r0$group=="rich",]$ego_id),NA,mean(node_r0$initial_degree),mean(node_r0$initial_wealth),gini(node_r0$initial_wealth),gmd(node_r0$initial_wealth),NA,mean(node_r0[node_r0$group=="rich",]$initial_degree),mean(node_r0[node_r0$group=="rich",]$initial_wealth),gini(node_r0[node_r0$group=="rich",]$initial_wealth),gmd(node_r0[node_r0$group=="rich",]$initial_wealth),NA,mean(node_r0[node_r0$group=="poor",]$initial_degree),mean(node_r0[node_r0$group=="poor",]$initial_wealth),gini(node_r0[node_r0$group=="poor",]$initial_wealth),gmd(node_r0[node_r0$group=="poor",]$initial_wealth),
                                       as.numeric(ifelse(is.na(table(node_r0$initial_degree<=isolationDegree)["TRUE"]),0,1)),
                                       as.numeric(sum(node_r0$everIsolated)/length(node_r0$ego_id)),
                                       NA,
                                       NA,
                                       NA,NA,NA
                                       )
      
      #For the loop at the next round (for round 1, the initial one is the same as the previous [1 prior] one)
      node_import = node_r0
      node_import$initial_coop = NA
      node_import$prev_coop = NA
      node_import$prev_wealth = node_import$initial_wealth
      node_import$prev_degree = node_import$initial_degree
      node_import$prev_avg_env_wealth = node_import$initial_avg_env_wealth
      node_import$prev_local_gini = node_import$initial_local_gini
      node_import$prev_rel_rank = node_import$initial_rel_rank
      node_import$prev_local_rate_coop = NA
      link_import = link_r0
      
      
      #####################################################
      # Section 3: Rounds 1 to 10 or more (behaviors in simulation: the equation of cooperation is different at round 1 because of no history)
      #3-1: Cooperation phase
      for (k in 1:L)
      {
        node_rX = node_import #Importing data
        node_rX$round = node_rX$round + 1
        node_rX[is.na(node_rX$prev_degree)==1,"prev_degree"] = 0
        node_rX[is.na(node_rX$prev_local_rate_coop)==1,"prev_local_rate_coop"] = 0
        #Only this calculation needs to change from Round 1
        if (k==1) {
          node_rX$prob_coop = as.numeric(V==0)*inv.logit((-1.816665) + (2.086067)*coop_rp1 + (1.800153)*coop_rp2) + as.numeric(V==1)*inv.logit((-2.031577) + (2.427157)*coop_rp1 + (1.684193)*coop_rp2 + (-1.528851)*GINI)
          } else {
          node_rX$prob_coop = as.numeric(V==0 & node_rX$prev_coop==0)*inv.logit(-1.039916) + as.numeric(V==0 & node_rX$prev_coop==1)*inv.logit(2.062023) + as.numeric(V==1 & node_rX$prev_coop==0)*inv.logit((-0.2574838)*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-1.214198)*GINI + (2.508148)*GINI*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-0.9749075)) + as.numeric(V==1 & node_rX$prev_coop==1)*inv.logit((- 0.6197254)*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-0.7480261)*GINI + (1.169674)*GINI*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (1.356784))
          } 
        #####manipulate cooperation rate ar round 1 depending on the degree! (keep the total coopertion rate at round 1 constant)
        #####make it so that the nodes with the top fractionCoop*100 percentile of degrees will 100% be a cooperator, but the average percentage of being a cooperator will be the same
        if(k==1){
          prob_coop_df = NULL
          nodesCoop = NULL
          nodesCoop = node_rX$prev_degree<=quantile(node_rX$prev_degree,fractionCoop)
          prob_coop_df = 
            data.frame(
              prob_coop = rev(node_rX$prob_coop[order(coop_rp)]),
              node_number = c(which(nodesCoop),which(!nodesCoop))
            )
          node_rX$prob_coop = prob_coop_df[order(prob_coop_df$node_number),]$prob_coop
          node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
        } else {
          node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
        }
        if (k==1) {
          node_rX$initial_coop = node_rX$coop
          } else {
          node_rX$initial_coop = node_rX$initial_coop
          }
        node_rX$cost = (-50)*node_rX$coop*node_rX$prev_degree
        node_rX$n_coop_received = NA
        for (i in 1:(dim(node_rX)[1]))
          {
          node_rX[i,]$n_coop_received = sum(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) ==
          "alt_id"]],"coop"])
          }
        node_rX$benefit = 100*node_rX$n_coop_received
        node_rX$payoff = node_rX$cost + node_rX$benefit
        node_rX$wealth = node_rX$prev_wealth + node_rX$payoff
        node_rX$rel_rank = NA
        node_rX$local_rate_coop = NA
        for (i in 1:dim(node_rX)[1])
          {
          node_rX[i,]$rel_rank = rank1(na.omit(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in%
          c("ego_id","alt_id")]],"wealth"]))/length(na.omit(node_rX[node_rX$ego_id %in%
          node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
          node_rX[i,]$local_rate_coop = mean(na.omit(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in%
          c("ego_id","alt_id")]],"coop"]))
          }
        node_rX$growth = as.numeric((node_rX$wealth/node_rX$prev_wealth) > 1)
        node_rX = node_rX[,c("ego_id","round","group","prev_degree","initial_wealth","initial_local_gini","initial_coop","coop","wealth","rel_rank","local_rate_coop","growth","everIsolated","maxDegreeLost")] #Pruning the previous-round data (degree is not updating yet)
        
        #3-2: Rewiring phase
        # 30% of ties (unidirectional) are being rewired
        link_rX_1 = link_import #Importing data (bidirectioanl ego-alter [ego_id < alter_id])
        colnames(link_rX_1) = c("ego_id","alt_id","prev_connected")
        link_rX_1$challenge = sample(0:1,dim(link_rX_1)[1],replace=TRUE,prob=c(0.7,0.3)) # The bidirectional ties being rewired are selected (rewiring rate = 0.3).
        ego_node_data =
        node_rX[,c("ego_id","wealth","coop","prev_degree","initial_wealth","initial_local_gini","initial_coop","rel_rank","local_rate_coop","growth")]
        colnames(ego_node_data) =
        c("ego_id","ego_wealth","ego_coop","ego_prev_degree","ego_initial_wealth","ego_initial_local_gini","ego_initial_coop","ego_rel_rank","ego_local_rate_coop","ego_growth")
        alt_node_data =
        node_rX[,c("ego_id","wealth","coop","prev_degree","initial_wealth","initial_local_gini","initial_coop","rel_rank","local_rate_coop","growth")]
        colnames(alt_node_data) =
        c("alt_id","alt_wealth","alt_coop","alt_prev_degree","alt_initial_wealth","alt_initial_local_gini","alt_initial_coop","alt_rel_rank","alt_local_rate_coop","alt_growth")
        link_rX_2 = merge(x=link_rX_1,y=ego_node_data,all.x=TRUE,all.y=FALSE,by="ego_id")
        link_rX_3 = merge(x=link_rX_2,y=alt_node_data,all.x=TRUE,all.y=FALSE,by="alt_id")
        link_rX_3$choice = sample(c("ego","alt"),dim(link_rX_3)[1],replace=TRUE,prob=c(0.5,0.5)) #decision maker for breaking a link, which is a unilateral decision
        #ego_prob: probability of choosing to connect when challenged (asked)
        link_rX_3$ego_prob = inv.logit((0.5134401)*link_rX_3$prev_connected + (-0.852406)*link_rX_3$ego_coop + (2.96549)*link_rX_3$alt_coop + (-0.1808545))
        link_rX_3$alt_prob = inv.logit((0.5134401)*link_rX_3$prev_connected + (-0.852406)*link_rX_3$alt_coop + (2.96549)*link_rX_3$ego_coop + (-0.1808545))
        link_rX_3$prob_connect = ifelse(link_rX_3$prev_connected == 1, ifelse(link_rX_3$choice == "ego", link_rX_3$ego_prob,
        link_rX_3$alt_prob), link_rX_3$ego_prob*link_rX_3$alt_prob)
        link_rX_3$connect_update = apply(data.frame(link_rX_3$prob_connect),1, function(x) {sample(1:0,1,prob=c(x,(1-x)))})
        link_rX_3$connected = ifelse(link_rX_3$challenge==0,link_rX_3$prev_connected,link_rX_3$connect_update)
        link_rX = link_rX_3[,c("ego_id","alt_id","connected")] #pruning and data is updated
        #Reflect the degree and local gini coefficient into the node data
        link_rXc_ego = link_rX[link_rX$connected==1,]
        link_rXc_alt = link_rX[link_rX$connected==1,]
        colnames(link_rXc_alt) = c("alt_id","ego_id","connected")
        link_rXc = rbind(link_rXc_ego,link_rXc_alt)
        link_rXc = link_rXc[order(link_rXc$ego_id),]
        link_rXc$alternumber = NA
        link_rXc[1,]$alternumber = 1
        for (i in 1:(dim(link_rXc)[1]-1))
          {
            if (link_rXc[i,]$ego_id == link_rXc[i+1,]$ego_id)
            {
            link_rXc[i+1,]$alternumber = link_rXc[i,]$alternumber + 1
            }
            else
            {
            link_rXc[i+1,]$alternumber = 1
            }
            #print(i)
          }
        link_rXc2 = reshape(link_rXc, direction = "wide", idvar=c("ego_id","connected"), timevar="alternumber")
        link_rXc2$degree = apply(link_rXc2[,colnames(link_rXc2)[substr(colnames(link_rXc2),1,3) == "alt"]],1,function(x) {length(na.omit(x))})
        node_rX_final = merge(x=node_rX[,c("ego_id","round","group","initial_wealth","initial_local_gini","initial_coop","coop","wealth","growth","everIsolated","maxDegreeLost")],y=link_rXc2,all.x=TRUE,all.y=FALSE,by="ego_id")
        node_rX_final[is.na(node_rX_final$degree)==1,"degree"] = 0
        node_rX_final$avg_env_wealth = NA
        node_rX_final$local_gini = NA #needs to be updated because the social network changes at the rewiring phase
        node_rX_final$local_rate_coop = NA
        node_rX_final$rel_rank = NA
        for (i in 1:dim(node_rX_final)[1])
          {
          node_rX_final[i,]$avg_env_wealth = mean(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
          node_rX_final[i,]$local_gini = gini(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
          node_rX_final[i,]$local_rate_coop = mean(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"coop"]))
          node_rX_final[i,]$rel_rank = rank1(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in%
          c("ego_id","alt_id")]],"wealth"]))/length(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
          node_rX_final[i,]$everIsolated = ifelse(node_rX_final[i,]$everIsolated==1,1,ifelse(node_rX_final[i,]$degree<=isolationDegree,1,0))
          node_rX_final[i,]$maxDegreeLost = pmax(node_r0[i,]$initial_degree - node_rX_final[i,]$degree, node_rX_final[i,]$maxDegreeLost, na.rm=TRUE)
        }
        
        
        #Finalization of round X and Visualization
        #plot(graph.data.frame(link_rX[link_rX$connected==1,],directed=F)) #plot.igraph
        result[result$round==k,2:25] =
        c(length(node_rX_final$ego_id),length(node_rX_final[node_rX_final$group=="rich",]$ego_id),mean(node_rX_final$coop),mean(node_rX_final$degree),mean(node_rX_final$wealth),gini(node_rX_final$wealth),gmd(node_rX_final$wealth),mean(node_rX_final[node_rX_final$group=="rich",]$coop),mean(node_rX_final[node_rX_final$group=="rich",]$degree),mean(node_rX_final[node_rX_final$group=="rich",]$wealth),gini(node_rX_final[node_rX_final$group=="rich",]$wealth),gmd(node_rX_final[node_rX_final$group=="rich",]$wealth),mean(node_rX_final[node_rX_final$group=="poor",]$coop),mean(node_rX_final[node_rX_final$group=="poor",]$degree),mean(node_rX_final[node_rX_final$group=="poor",]$wealth),gini(node_rX_final[node_rX_final$group=="poor",]$wealth),gmd(node_rX_final[node_rX_final$group=="poor",]$wealth),
                                       as.numeric(ifelse(is.na(table(node_rX_final$degree<=isolationDegree)["TRUE"]),0,1)),
                                       as.numeric(sum(node_rX_final$everIsolated)/length(node_rX_final$ego_id)),
          prop.table(table(node_rX_final[node_rX_final$initial_coop==1]$coop))["0"],
          prop.table(table(node_rX_final[node_rX_final$initial_coop==0]$coop))["1"],
          suppressWarnings({mean(node_rX_final$maxDegreeLost,na.rm=TRUE)}),
          suppressWarnings({mean(node_rX_final[node_rX_final$initial_coop==1]$maxDegreeLost,na.rm=TRUE)}),
          suppressWarnings({mean(node_rX_final[node_rX_final$initial_coop==0]$maxDegreeLost,na.rm=TRUE)})
          )
        
        #For the loop
        node_import = node_rX_final
        colnames(node_import)[colnames(node_import) %in%
        c("coop","wealth","growth","degree","avg_env_wealth","local_gini","local_rate_coop","rel_rank")] =
        c("prev_coop","prev_wealth","prev_growth","prev_degree","prev_avg_env_wealth","prev_local_gini","prev_local_rate_coop","prev_rel_rank")
        link_import = link_rX
        #print(paste0("Round ",k," is done."))
      }
      
      coopFrac[m] = fractionCoop
      avgCoop[m] = result[result$round==1,]$avg_coop
      avgCoopFinal[m] = result[result$round==10,]$avg_coop
      percentIsolation[m] = max(result[result$round>=1,]$percentIsolation)
      isolation[m] = max(result[result$round>=1,]$isolation)
      percentIsolatedD[m] = prop.table(table(node_rX_final[node_rX_final$everIsolated==1,]$initial_coop))["0"]
      
      
      link_rX_final = data.table::melt(setDT(node_rX_final),
                               measure = patterns('alt_id'),
                               variable.name = 'linkNumber', 
                               value.name = c('alt_id'))
      link_rX_final = data.frame(link_rX_final)[c("ego_id","alt_id")]
      link_rX_final = link_rX_final[complete.cases(link_rX_final),]
      link_rX_final = data.frame(t(unique(apply(link_rX_final, 1, function(x) sort(x))))) %>% distinct(X1, X2)
      node_g_final = data.frame(node_rX_final)[c("ego_id","initial_coop","coop")]
      node_g_final$initial_coop = factor(node_g_final$initial_coop)
      g_rX_final = graph_from_data_frame(link_rX_final, directed = FALSE, vertices=node_g_final)
      g_r0 = graph_from_data_frame(link_r0[link_r0$connected==1,][1:2], directed = FALSE, vertices=node_r0)
      nCommunities[m] = max(membership(cluster_louvain(g_rX_final)),na.rm=TRUE)
      communitySize[m] = mean(table(membership(cluster_louvain(g_rX_final))),na.rm=TRUE)
      assortativityInitial[m] = assortativity(g_r0, V(g_rX_final)$initial_coop == 1)
      assortativityFinal[m] = assortativity(g_rX_final, V(g_rX_final)$initial_coop == 1)
      conversionRate[m] = prop.table(table(V(g_rX_final)$coop == V(g_rX_final)$initial_coop))["FALSE"]
      conversionToD[m] = prop.table(table(V(g_rX_final)$coop[V(g_rX_final)$initial_coop==1]))["0"]
      conversionToC[m] = prop.table(table(V(g_rX_final)$coop[V(g_rX_final)$initial_coop==0]))["1"]
      degreeC[m] = mean(igraph::degree(g_rX_final)[V(g_rX_final)$coop==1],na.rm=TRUE)
      degreeD[m] = mean(igraph::degree(g_rX_final)[V(g_rX_final)$coop==0],na.rm=TRUE)
      meanConversionToD[m] = mean(result[result$round>=2,]$meanConversionToD, na.rm=TRUE)
      meanConversionToC[m] = mean(result[result$round>=2,]$meanConversionToC, na.rm=TRUE)
      degreeLost[m] = result[result$round==10,]$degreeLost
      degreeLostC[m] = result[result$round==10,]$degreeLostC
      degreeLostD[m] = result[result$round==10,]$degreeLostD
  
    }
    
df.netIntLowDegree = rbind(df.netIntLowDegree,
                  data.frame(
                    coopFrac = coopFrac,
                    avgCoop = avgCoop,
                    avgCoopFinal = avgCoopFinal,
                    percentIsolation = percentIsolation,
                    isolation = isolation,
                    percentIsolatedD = percentIsolatedD,
                    nCommunities = nCommunities,
                    communitySize = communitySize,
                    assortativityInitial = assortativityInitial,
                    assortativityFinal = assortativityFinal,
                    conversionRate = conversionRate,
                    conversionToD = conversionToD, 
                    conversionToC = conversionToC, 
                    degreeC = degreeC, 
                    degreeD = degreeD,
                    meanConversionToD = meanConversionToD, 
                    meanConversionToC = meanConversionToC,
                    degreeLost = degreeLost,
                    degreeLostC = degreeLostC,
                    degreeLostD = degreeLostD
                    ))
plot(g_r0,vertex.color=V(g_rX_final)$initial_coop,vertex.label=ifelse(is.na(V(g_rX_final)$initial_coop),"NA",ifelse(V(g_rX_final)$initial_coop==1,"C","D")),main=paste("fracCoop=",frac,", round 0",sep=""))
plot(g_rX_final,vertex.color=V(g_rX_final)$initial_coop,vertex.label=ifelse(is.na(V(g_rX_final)$initial_coop),"NA",ifelse(V(g_rX_final)$initial_coop==1,"C","D")),main=paste("fracCoop=",frac,", final round",sep=""))
}

sum.netIntLowDegree <- data.frame(
  df.netIntLowDegree %>% 
    group_by(coopFrac) %>%
      summarise(
        mean.isolation = mean(isolation),
        ci.isolation   = 1.96 * sd(isolation)/sqrt(n()),
        mean.percentIsolation = mean(percentIsolation),
        ci.percentIsolation   = 1.96 * sd(percentIsolation)/sqrt(n()),
        mean.percentIsolatedD = mean(percentIsolatedD,na.rm=TRUE),
        ci.percentIsolatedD   = 1.96 * sd(percentIsolatedD,na.rm=TRUE)/sqrt(sum(isolation)),
        mean.avgCoop = mean(avgCoop,na.rm=TRUE),
        ci.avgCoop   = 1.96 * sd(avgCoop,na.rm=TRUE)/sqrt(n()),
        mean.avgCoopFinal = mean(avgCoopFinal,na.rm=TRUE),
        ci.avgCoopFinal   = 1.96 * sd(avgCoopFinal,na.rm=TRUE)/sqrt(n()),
        mean.nCommunities = mean(nCommunities,na.rm=TRUE),
        ci.nCommunities   = 1.96 * sd(nCommunities,na.rm=TRUE)/sqrt(n()),
        mean.communitySize = mean(communitySize,na.rm=TRUE),
        ci.communitySize   = 1.96 * sd(communitySize,na.rm=TRUE)/sqrt(n()),
        mean.assortativityInitial = mean(assortativityInitial,na.rm=TRUE),
        ci.assortativityInitial   = 1.96 * sd(assortativityInitial,na.rm=TRUE)/sqrt(n()),
        mean.assortativityFinal = mean(assortativityFinal,na.rm=TRUE),
        ci.assortativityFinal   = 1.96 * sd(assortativityFinal,na.rm=TRUE)/sqrt(n()),
        mean.conversionRate = mean(conversionRate,na.rm=TRUE),
        ci.conversionRate   = 1.96 * sd(conversionRate,na.rm=TRUE)/sqrt(n()),
        mean.conversionToD = mean(conversionToD,na.rm=TRUE),
        ci.conversionToD   = 1.96 * sd(conversionToD,na.rm=TRUE)/sqrt(n()),
        mean.conversionToC = mean(conversionToC,na.rm=TRUE),
        ci.conversionToC   = 1.96 * sd(conversionToC,na.rm=TRUE)/sqrt(n()),
        mean.degreeC = mean(degreeC,na.rm=TRUE),
        ci.degreeC   = 1.96 * sd(degreeC,na.rm=TRUE)/sqrt(n()),
        mean.degreeD = mean(degreeD,na.rm=TRUE),
        ci.degreeD   = 1.96 * sd(degreeD,na.rm=TRUE)/sqrt(n()),
        mean.meanConversionToD = mean(meanConversionToD,na.rm=TRUE),
        ci.meanConversionToD   = 1.96 * sd(meanConversionToD,na.rm=TRUE)/sqrt(n()),
        mean.meanConversionToC = mean(meanConversionToC,na.rm=TRUE),
        ci.meanConversionToC   = 1.96 * sd(meanConversionToC,na.rm=TRUE)/sqrt(n()),
        mean.degreeLost = mean(degreeLost,na.rm=TRUE),
        ci.degreeLost   = 1.96 * sd(degreeLost,na.rm=TRUE)/sqrt(n()),
        mean.degreeLostC = mean(degreeLostC,na.rm=TRUE),
        ci.degreeLostC   = 1.96 * sd(degreeLostC,na.rm=TRUE)/sqrt(n()),
        mean.degreeLostD = mean(degreeLostD,na.rm=TRUE),
        ci.degreeLostD   = 1.96 * sd(degreeLostD,na.rm=TRUE)/sqrt(n())
        )
  )

kable(sum.netIntLowDegree[c(1:9)]) %>% kableExtra::kable_styling(font_size = 10)

coopFrac	mean.isolation	ci.isolation	mean.percentIsolation	ci.percentIsolation	mean.percentIsolatedD	ci.percentIsolatedD	mean.avgCoop	ci.avgCoop
0.0	0.532	0.0437809	0.0460000	0.0046922	0.8752159	0.0259206	0.6969412	0.0099184
0.1	0.486	0.0438536	0.0410588	0.0044243	0.8644689	0.0280003	0.6992941	0.0097812
0.2	0.468	0.0437809	0.0414118	0.0048647	0.8301075	0.0317277	0.7015294	0.0102715
0.3	0.434	0.0434869	0.0375294	0.0044964	0.8370370	0.0312692	0.7016471	0.0096155
0.4	0.430	0.0434388	0.0349412	0.0042597	0.8855634	0.0280366	0.7040000	0.0099747
0.5	0.432	0.0434632	0.0345882	0.0041104	0.8795662	0.0299354	0.7032941	0.0094693
0.6	0.476	0.0438203	0.0387059	0.0042593	0.8581224	0.0294044	0.6997647	0.0100892
0.7	0.486	0.0438536	0.0418824	0.0045026	0.8554913	0.0290239	0.6934118	0.0098512
0.8	0.490	0.0438621	0.0407059	0.0044193	0.8905678	0.0255032	0.6984706	0.0096146
0.9	0.518	0.0438424	0.0467059	0.0049922	0.8835846	0.0253330	0.6921176	0.0100919
1.0	0.556	0.0435948	0.0529412	0.0054962	0.8781061	0.0247481	0.6911765	0.0102598

kable(sum.netIntLowDegree[c(1,10:17)]) %>% kableExtra::kable_styling(font_size = 10) 

coopFrac	mean.avgCoopFinal	ci.avgCoopFinal	mean.nCommunities	ci.nCommunities	mean.communitySize	ci.communitySize	mean.assortativityInitial	ci.assortativityInitial
0.0	0.6924706	0.0097983	2.682	0.0468310	6.615833	0.1238711	-0.0704514	0.0110977
0.1	0.7023529	0.0101303	2.678	0.0456186	6.615833	0.1223652	-0.0701275	0.0111441
0.2	0.7005882	0.0100159	2.706	0.0447021	6.536500	0.1198070	-0.0721554	0.0111999
0.3	0.6883529	0.0099188	2.666	0.0452907	6.644167	0.1224157	-0.0748324	0.0108500
0.4	0.6885882	0.0100235	2.658	0.0458475	6.669667	0.1234736	-0.0790560	0.0099904
0.5	0.6972941	0.0100196	2.676	0.0453399	6.618667	0.1220691	-0.0844428	0.0110704
0.6	0.6960000	0.0100291	2.700	0.0445672	6.550667	0.1199114	-0.0800918	0.0103097
0.7	0.7024706	0.0096374	2.638	0.0449932	6.715000	0.1230765	-0.0878133	0.0110702
0.8	0.6978824	0.0097349	2.688	0.0466482	6.598833	0.1233671	-0.0718445	0.0110094
0.9	0.7022353	0.0095143	2.708	0.0432288	6.519500	0.1176331	-0.0726472	0.0109855
1.0	0.6940000	0.0095465	2.690	0.0452444	6.581833	0.1213243	-0.0753249	0.0112991

kable(sum.netIntLowDegree[c(1,18:25)]) %>% kableExtra::kable_styling(font_size = 10) 

coopFrac	mean.assortativityFinal	ci.assortativityFinal	mean.conversionRate	ci.conversionRate	mean.conversionToD	ci.conversionToD	mean.conversionToC	ci.conversionToC
0.0	-0.0588831	0.0054180	0.4171765	0.0102754	0.3079011	0.0112816	0.6938654	0.0180888
0.1	-0.0550204	0.0054722	0.4200000	0.0107081	0.3079973	0.0115326	0.7213619	0.0181675
0.2	-0.0627492	0.0052061	0.4047059	0.0100792	0.2923497	0.0109386	0.6992634	0.0183045
0.3	-0.0578321	0.0053261	0.4189412	0.0101811	0.3142649	0.0113727	0.6904694	0.0185433
0.4	-0.0556872	0.0053338	0.4095294	0.0103233	0.3119411	0.0109740	0.6868707	0.0190884
0.5	-0.0601939	0.0053717	0.4194118	0.0103497	0.3101919	0.0111733	0.7030567	0.0186047
0.6	-0.0572647	0.0054577	0.4171765	0.0104042	0.3078311	0.0111042	0.6992666	0.0188132
0.7	-0.0561214	0.0051075	0.4189412	0.0104446	0.3033502	0.0114476	0.7055446	0.0179022
0.8	-0.0562569	0.0054781	0.4109412	0.0106750	0.3022745	0.0111135	0.6868295	0.0181381
0.9	-0.0609449	0.0055895	0.4216471	0.0100715	0.3008787	0.0109256	0.7137053	0.0178200
1.0	-0.0596674	0.0055195	0.4138824	0.0101123	0.3015647	0.0111019	0.6865431	0.0175059

kable(sum.netIntLowDegree[c(1,26:33)]) %>% kableExtra::kable_styling(font_size = 10) 

coopFrac	mean.degreeC	ci.degreeC	mean.degreeD	ci.degreeD	mean.meanConversionToD	ci.meanConversionToD	mean.meanConversionToC	ci.meanConversionToC
0.0	11.02560	0.0843788	8.728804	0.1168585	0.3023525	0.0060099	0.6789940	0.0061272
0.1	11.11527	0.0842879	8.745885	0.1198110	0.3018477	0.0058163	0.6950404	0.0064098
0.2	11.08739	0.0807612	8.699833	0.1178646	0.3030937	0.0059827	0.6862745	0.0064254
0.3	11.10344	0.0870042	8.780629	0.1151682	0.2992981	0.0059488	0.6738834	0.0068023
0.4	11.16911	0.0888452	8.755639	0.1180366	0.2949554	0.0063474	0.6835372	0.0058913
0.5	11.11939	0.0801389	8.739078	0.1167014	0.3001828	0.0060898	0.6910638	0.0056423
0.6	11.11432	0.0831961	8.674462	0.1176464	0.2995744	0.0056493	0.6865402	0.0060294
0.7	11.10209	0.0869425	8.691486	0.1211137	0.3014602	0.0060422	0.6944444	0.0059223
0.8	11.08310	0.0835440	8.727968	0.1177927	0.2959288	0.0056759	0.6788913	0.0060580
0.9	11.06081	0.0772334	8.690480	0.1152031	0.3072473	0.0058074	0.6898005	0.0064664
1.0	11.06020	0.0842639	8.670291	0.1158618	0.3031632	0.0061387	0.6806997	0.0064272

kable(sum.netIntLowDegree[c(1,34:ncol(sum.netIntLowDegree))]) %>% kableExtra::kable_styling(font_size = 10) 

coopFrac	mean.degreeLost	ci.degreeLost	mean.degreeLostC	ci.degreeLostC	mean.degreeLostD	ci.degreeLostD
0.0	-0.8638554	0.0547047	-0.9673806	0.0541616	-0.6592700	0.0513858
0.1	-0.8815420	0.0545608	-0.9805756	0.0548255	-0.7129173	0.0509411
0.2	-0.8487471	0.0546166	-0.9214822	0.0541931	-0.6898417	0.0530895
0.3	-0.8815809	0.0532822	-0.9589127	0.0528883	-0.7314663	0.0517045
0.4	-0.8761838	0.0575615	-0.9176898	0.0591551	-0.7977304	0.0539197
0.5	-0.9512353	0.0571366	-1.0118534	0.0591979	-0.8124516	0.0502109
0.6	-0.8443456	0.0529607	-0.8880107	0.0519384	-0.7564892	0.0543042
0.7	-0.8731397	0.0531126	-0.9187655	0.0533836	-0.7614351	0.0513190
0.8	-0.8635882	0.0547863	-0.9366237	0.0554279	-0.6661220	0.0491773
0.9	-0.8469632	0.0540782	-0.9520010	0.0519941	-0.6216935	0.0533232
1.0	-0.8030956	0.0511954	-0.8897691	0.0492045	-0.6027707	0.0515510

compare_means(percentIsolation ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.              group1 group2        p p.adj p.format p.signif method  
##    <chr>            <chr>  <chr>     <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 percentIsolation 0      0.1    0.135    1     0.13477  ns       Wilcoxon
##  2 percentIsolation 0      0.2    0.0515   1     0.05153  ns       Wilcoxon
##  3 percentIsolation 0      0.3    0.00307  0.14  0.00307  **       Wilcoxon
##  4 percentIsolation 0      0.4    0.000324 0.017 0.00032  ***      Wilcoxon
##  5 percentIsolation 0      0.5    0.000391 0.02  0.00039  ***      Wilcoxon
##  6 percentIsolation 0      0.6    0.0322   1     0.03225  *        Wilcoxon
##  7 percentIsolation 0      0.7    0.187    1     0.18674  ns       Wilcoxon
##  8 percentIsolation 0      0.8    0.108    1     0.10810  ns       Wilcoxon
##  9 percentIsolation 0      0.9    0.841    1     0.84132  ns       Wilcoxon
## 10 percentIsolation 0      1      0.247    1     0.24663  ns       Wilcoxon
## # … with 45 more rows

compare_means(avgCoop ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.     group1 group2     p p.adj p.format p.signif method  
##    <chr>   <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 avgCoop 0      0.1    0.673     1 0.673    ns       Wilcoxon
##  2 avgCoop 0      0.2    0.486     1 0.486    ns       Wilcoxon
##  3 avgCoop 0      0.3    0.498     1 0.498    ns       Wilcoxon
##  4 avgCoop 0      0.4    0.315     1 0.315    ns       Wilcoxon
##  5 avgCoop 0      0.5    0.393     1 0.393    ns       Wilcoxon
##  6 avgCoop 0      0.6    0.663     1 0.663    ns       Wilcoxon
##  7 avgCoop 0      0.7    0.736     1 0.736    ns       Wilcoxon
##  8 avgCoop 0      0.8    0.648     1 0.648    ns       Wilcoxon
##  9 avgCoop 0      0.9    0.420     1 0.420    ns       Wilcoxon
## 10 avgCoop 0      1      0.435     1 0.435    ns       Wilcoxon
## # … with 45 more rows

compare_means(avgCoopFinal ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.          group1 group2     p p.adj p.format p.signif method  
##    <chr>        <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 avgCoopFinal 0      0.1    0.149     1 0.149    ns       Wilcoxon
##  2 avgCoopFinal 0      0.2    0.195     1 0.195    ns       Wilcoxon
##  3 avgCoopFinal 0      0.3    0.629     1 0.629    ns       Wilcoxon
##  4 avgCoopFinal 0      0.4    0.536     1 0.536    ns       Wilcoxon
##  5 avgCoopFinal 0      0.5    0.427     1 0.427    ns       Wilcoxon
##  6 avgCoopFinal 0      0.6    0.619     1 0.619    ns       Wilcoxon
##  7 avgCoopFinal 0      0.7    0.180     1 0.180    ns       Wilcoxon
##  8 avgCoopFinal 0      0.8    0.532     1 0.532    ns       Wilcoxon
##  9 avgCoopFinal 0      0.9    0.164     1 0.164    ns       Wilcoxon
## 10 avgCoopFinal 0      1      0.762     1 0.762    ns       Wilcoxon
## # … with 45 more rows

compare_means(nCommunities ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.          group1 group2     p p.adj p.format p.signif method  
##    <chr>        <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 nCommunities 0      0.1    0.966     1 0.966    ns       Wilcoxon
##  2 nCommunities 0      0.2    0.402     1 0.402    ns       Wilcoxon
##  3 nCommunities 0      0.3    0.707     1 0.707    ns       Wilcoxon
##  4 nCommunities 0      0.4    0.520     1 0.520    ns       Wilcoxon
##  5 nCommunities 0      0.5    0.933     1 0.933    ns       Wilcoxon
##  6 nCommunities 0      0.6    0.503     1 0.503    ns       Wilcoxon
##  7 nCommunities 0      0.7    0.237     1 0.237    ns       Wilcoxon
##  8 nCommunities 0      0.8    0.852     1 0.852    ns       Wilcoxon
##  9 nCommunities 0      0.9    0.322     1 0.322    ns       Wilcoxon
## 10 nCommunities 0      1      0.739     1 0.739    ns       Wilcoxon
## # … with 45 more rows

compare_means(communitySize ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.           group1 group2     p p.adj p.format p.signif method  
##    <chr>         <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 communitySize 0      0.1    0.966     1 0.966    ns       Wilcoxon
##  2 communitySize 0      0.2    0.402     1 0.402    ns       Wilcoxon
##  3 communitySize 0      0.3    0.707     1 0.707    ns       Wilcoxon
##  4 communitySize 0      0.4    0.520     1 0.520    ns       Wilcoxon
##  5 communitySize 0      0.5    0.933     1 0.933    ns       Wilcoxon
##  6 communitySize 0      0.6    0.503     1 0.503    ns       Wilcoxon
##  7 communitySize 0      0.7    0.237     1 0.237    ns       Wilcoxon
##  8 communitySize 0      0.8    0.852     1 0.852    ns       Wilcoxon
##  9 communitySize 0      0.9    0.322     1 0.322    ns       Wilcoxon
## 10 communitySize 0      1      0.739     1 0.739    ns       Wilcoxon
## # … with 45 more rows

compare_means(assortativityInitial ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.                  group1 group2     p p.adj p.format p.signif method  
##    <chr>                <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 assortativityInitial 0      0.1    0.720     1 0.720    ns       Wilcoxon
##  2 assortativityInitial 0      0.2    0.845     1 0.845    ns       Wilcoxon
##  3 assortativityInitial 0      0.3    0.746     1 0.746    ns       Wilcoxon
##  4 assortativityInitial 0      0.4    0.416     1 0.416    ns       Wilcoxon
##  5 assortativityInitial 0      0.5    0.110     1 0.110    ns       Wilcoxon
##  6 assortativityInitial 0      0.6    0.376     1 0.376    ns       Wilcoxon
##  7 assortativityInitial 0      0.7    0.106     1 0.106    ns       Wilcoxon
##  8 assortativityInitial 0      0.8    0.874     1 0.874    ns       Wilcoxon
##  9 assortativityInitial 0      0.9    0.924     1 0.924    ns       Wilcoxon
## 10 assortativityInitial 0      1      0.640     1 0.640    ns       Wilcoxon
## # … with 45 more rows

compare_means(assortativityFinal ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.                group1 group2     p p.adj p.format p.signif method  
##    <chr>              <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 assortativityFinal 0      0.1    0.467     1 0.467    ns       Wilcoxon
##  2 assortativityFinal 0      0.2    0.238     1 0.238    ns       Wilcoxon
##  3 assortativityFinal 0      0.3    0.791     1 0.791    ns       Wilcoxon
##  4 assortativityFinal 0      0.4    0.355     1 0.355    ns       Wilcoxon
##  5 assortativityFinal 0      0.5    0.721     1 0.721    ns       Wilcoxon
##  6 assortativityFinal 0      0.6    0.681     1 0.681    ns       Wilcoxon
##  7 assortativityFinal 0      0.7    0.471     1 0.471    ns       Wilcoxon
##  8 assortativityFinal 0      0.8    0.711     1 0.711    ns       Wilcoxon
##  9 assortativityFinal 0      0.9    0.582     1 0.582    ns       Wilcoxon
## 10 assortativityFinal 0      1      0.679     1 0.679    ns       Wilcoxon
## # … with 45 more rows

compare_means(conversionRate ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.            group1 group2     p p.adj p.format p.signif method  
##    <chr>          <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 conversionRate 0      0.1    0.761     1 0.761    ns       Wilcoxon
##  2 conversionRate 0      0.2    0.116     1 0.116    ns       Wilcoxon
##  3 conversionRate 0      0.3    0.984     1 0.984    ns       Wilcoxon
##  4 conversionRate 0      0.4    0.336     1 0.336    ns       Wilcoxon
##  5 conversionRate 0      0.5    0.821     1 0.821    ns       Wilcoxon
##  6 conversionRate 0      0.6    0.952     1 0.952    ns       Wilcoxon
##  7 conversionRate 0      0.7    0.995     1 0.995    ns       Wilcoxon
##  8 conversionRate 0      0.8    0.310     1 0.310    ns       Wilcoxon
##  9 conversionRate 0      0.9    0.675     1 0.675    ns       Wilcoxon
## 10 conversionRate 0      1      0.778     1 0.778    ns       Wilcoxon
## # … with 45 more rows

compare_means(conversionToD ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.           group1 group2      p p.adj p.format p.signif method  
##    <chr>         <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 conversionToD 0      0.1    0.829      1 0.8286   ns       Wilcoxon
##  2 conversionToD 0      0.2    0.0511     1 0.0511   ns       Wilcoxon
##  3 conversionToD 0      0.3    0.562      1 0.5617   ns       Wilcoxon
##  4 conversionToD 0      0.4    0.587      1 0.5871   ns       Wilcoxon
##  5 conversionToD 0      0.5    0.863      1 0.8630   ns       Wilcoxon
##  6 conversionToD 0      0.6    0.924      1 0.9236   ns       Wilcoxon
##  7 conversionToD 0      0.7    0.403      1 0.4031   ns       Wilcoxon
##  8 conversionToD 0      0.8    0.514      1 0.5144   ns       Wilcoxon
##  9 conversionToD 0      0.9    0.399      1 0.3987   ns       Wilcoxon
## 10 conversionToD 0      1      0.400      1 0.3996   ns       Wilcoxon
## # … with 45 more rows

compare_means(conversionToC ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.           group1 group2      p p.adj p.format p.signif method  
##    <chr>         <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 conversionToC 0      0.1    0.0601     1 0.060    ns       Wilcoxon
##  2 conversionToC 0      0.2    0.754      1 0.754    ns       Wilcoxon
##  3 conversionToC 0      0.3    0.700      1 0.700    ns       Wilcoxon
##  4 conversionToC 0      0.4    0.489      1 0.489    ns       Wilcoxon
##  5 conversionToC 0      0.5    0.487      1 0.487    ns       Wilcoxon
##  6 conversionToC 0      0.6    0.708      1 0.708    ns       Wilcoxon
##  7 conversionToC 0      0.7    0.450      1 0.450    ns       Wilcoxon
##  8 conversionToC 0      0.8    0.522      1 0.522    ns       Wilcoxon
##  9 conversionToC 0      0.9    0.125      1 0.125    ns       Wilcoxon
## 10 conversionToC 0      1      0.481      1 0.481    ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeC ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.     group1 group2      p p.adj p.format p.signif method  
##    <chr>   <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeC 0      0.1    0.203   1    0.203    ns       Wilcoxon
##  2 degreeC 0      0.2    0.357   1    0.357    ns       Wilcoxon
##  3 degreeC 0      0.3    0.254   1    0.254    ns       Wilcoxon
##  4 degreeC 0      0.4    0.0171  0.94 0.017    *        Wilcoxon
##  5 degreeC 0      0.5    0.331   1    0.331    ns       Wilcoxon
##  6 degreeC 0      0.6    0.133   1    0.133    ns       Wilcoxon
##  7 degreeC 0      0.7    0.178   1    0.178    ns       Wilcoxon
##  8 degreeC 0      0.8    0.382   1    0.382    ns       Wilcoxon
##  9 degreeC 0      0.9    0.743   1    0.743    ns       Wilcoxon
## 10 degreeC 0      1      0.607   1    0.607    ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeD ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.     group1 group2     p p.adj p.format p.signif method  
##    <chr>   <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeD 0      0.1    0.696     1 0.70     ns       Wilcoxon
##  2 degreeD 0      0.2    0.534     1 0.53     ns       Wilcoxon
##  3 degreeD 0      0.3    0.513     1 0.51     ns       Wilcoxon
##  4 degreeD 0      0.4    0.718     1 0.72     ns       Wilcoxon
##  5 degreeD 0      0.5    0.938     1 0.94     ns       Wilcoxon
##  6 degreeD 0      0.6    0.460     1 0.46     ns       Wilcoxon
##  7 degreeD 0      0.7    0.849     1 0.85     ns       Wilcoxon
##  8 degreeD 0      0.8    0.787     1 0.79     ns       Wilcoxon
##  9 degreeD 0      0.9    0.542     1 0.54     ns       Wilcoxon
## 10 degreeD 0      1      0.509     1 0.51     ns       Wilcoxon
## # … with 45 more rows

compare_means(meanConversionToD ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.               group1 group2     p p.adj p.format p.signif method  
##    <chr>             <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 meanConversionToD 0      0.1    0.910     1 0.910    ns       Wilcoxon
##  2 meanConversionToD 0      0.2    0.939     1 0.939    ns       Wilcoxon
##  3 meanConversionToD 0      0.3    0.549     1 0.549    ns       Wilcoxon
##  4 meanConversionToD 0      0.4    0.155     1 0.155    ns       Wilcoxon
##  5 meanConversionToD 0      0.5    0.877     1 0.877    ns       Wilcoxon
##  6 meanConversionToD 0      0.6    0.690     1 0.690    ns       Wilcoxon
##  7 meanConversionToD 0      0.7    0.883     1 0.883    ns       Wilcoxon
##  8 meanConversionToD 0      0.8    0.230     1 0.230    ns       Wilcoxon
##  9 meanConversionToD 0      0.9    0.341     1 0.341    ns       Wilcoxon
## 10 meanConversionToD 0      1      0.997     1 0.997    ns       Wilcoxon
## # … with 45 more rows

compare_means(meanConversionToC ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.               group1 group2     p p.adj p.format p.signif method  
##    <chr>             <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 meanConversionToC 0      0.1    0.154     1 0.154    ns       Wilcoxon
##  2 meanConversionToC 0      0.2    0.358     1 0.358    ns       Wilcoxon
##  3 meanConversionToC 0      0.3    0.667     1 0.667    ns       Wilcoxon
##  4 meanConversionToC 0      0.4    0.680     1 0.680    ns       Wilcoxon
##  5 meanConversionToC 0      0.5    0.333     1 0.333    ns       Wilcoxon
##  6 meanConversionToC 0      0.6    0.357     1 0.357    ns       Wilcoxon
##  7 meanConversionToC 0      0.7    0.122     1 0.122    ns       Wilcoxon
##  8 meanConversionToC 0      0.8    0.892     1 0.892    ns       Wilcoxon
##  9 meanConversionToC 0      0.9    0.351     1 0.351    ns       Wilcoxon
## 10 meanConversionToC 0      1      0.837     1 0.837    ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeLost ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.        group1 group2      p p.adj p.format p.signif method  
##    <chr>      <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeLost 0      0.1    0.889      1 0.88876  ns       Wilcoxon
##  2 degreeLost 0      0.2    0.564      1 0.56378  ns       Wilcoxon
##  3 degreeLost 0      0.3    0.982      1 0.98217  ns       Wilcoxon
##  4 degreeLost 0      0.4    0.958      1 0.95833  ns       Wilcoxon
##  5 degreeLost 0      0.5    0.0684     1 0.06844  ns       Wilcoxon
##  6 degreeLost 0      0.6    0.573      1 0.57253  ns       Wilcoxon
##  7 degreeLost 0      0.7    0.890      1 0.89031  ns       Wilcoxon
##  8 degreeLost 0      0.8    0.839      1 0.83899  ns       Wilcoxon
##  9 degreeLost 0      0.9    0.457      1 0.45731  ns       Wilcoxon
## 10 degreeLost 0      1      0.0673     1 0.06726  ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeLostC ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.         group1 group2      p p.adj p.format p.signif method  
##    <chr>       <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeLostC 0      0.1    0.963      1 0.963    ns       Wilcoxon
##  2 degreeLostC 0      0.2    0.299      1 0.299    ns       Wilcoxon
##  3 degreeLostC 0      0.3    0.617      1 0.617    ns       Wilcoxon
##  4 degreeLostC 0      0.4    0.254      1 0.254    ns       Wilcoxon
##  5 degreeLostC 0      0.5    0.473      1 0.473    ns       Wilcoxon
##  6 degreeLostC 0      0.6    0.0930     1 0.093    ns       Wilcoxon
##  7 degreeLostC 0      0.7    0.327      1 0.327    ns       Wilcoxon
##  8 degreeLostC 0      0.8    0.511      1 0.511    ns       Wilcoxon
##  9 degreeLostC 0      0.9    0.529      1 0.529    ns       Wilcoxon
## 10 degreeLostC 0      1      0.0757     1 0.076    ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeLostD ~ coopFrac, data=df.netIntLowDegree)

## # A tibble: 55 × 8
##    .y.         group1 group2      p p.adj p.format p.signif method  
##    <chr>       <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeLostD 0      0.1    0.522      1 0.52183  ns       Wilcoxon
##  2 degreeLostD 0      0.2    0.858      1 0.85836  ns       Wilcoxon
##  3 degreeLostD 0      0.3    0.454      1 0.45450  ns       Wilcoxon
##  4 degreeLostD 0      0.4    0.0910     1 0.09098  ns       Wilcoxon
##  5 degreeLostD 0      0.5    0.0394     1 0.03941  *        Wilcoxon
##  6 degreeLostD 0      0.6    0.235      1 0.23516  ns       Wilcoxon
##  7 degreeLostD 0      0.7    0.199      1 0.19871  ns       Wilcoxon
##  8 degreeLostD 0      0.8    0.869      1 0.86895  ns       Wilcoxon
##  9 degreeLostD 0      0.9    0.418      1 0.41849  ns       Wilcoxon
## 10 degreeLostD 0      1      0.182      1 0.18246  ns       Wilcoxon
## # … with 45 more rows

g.netIntLowDegree = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="percentIsolation", add = "mean_se") +
  stat_compare_means(ref.group = "0", label = "p.signif", label.y = 0.098, method="t.test") +  
  labs(
    title = "Isolation when cooperators are assigned to low-degree nodes",
    x = "Proportion of low-degree nodes assigned to cooperators",
    y = "Propoption of ever-isolated individuals") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(ylim=c(0,0.10))
g.netIntLowDegree

Several explanations

More beneficial to break links

Mathematically, the expected benefit of breaking/not making new links only depends on the total number of cooperators and defectors, not the number of C and D around an ego

Altruism: making/retaining links even when there is no benefit

Degree changes in initial defectors are no different across conditions

Evolution of cooperation

This seems to occur regardless of initial defector/cooperator position (the same number of C's become D's, and the same number of D's become C's)

D has more links to start

This seems to be the most plausible explanation. Regardless of the network position, initial D's ends up becoming C's in most cases, but not without a cost: many D's tend to lose links first before reaching equilibrium, so networks that have the greatest D degrees may have the least isolation.

Simulation 2

All graphs start as Erdos-Renyi random newtorks

Change where to place the cooperators in the network: 1. cooperators have high degrees

#Define degrees of isolation
isolationDegree = 2

#number of iterations per arm
iterations = 500


df.netIntHighDegree = data.frame(
  coopFrac = NULL,
  avgCoop = NULL, 
  avgCoopFinal = NULL, 
  percentIsolation = NULL,
  isolation = NULL,
  percentIsolatedD = NULL,
  nCommunities = NULL,
  communitySize = NULL,
  assortativityInitial = NULL,
  assortativityFinal = NULL,
  conversionRate = NULL,
  conversionToD = NULL, 
  conversionToC = NULL, 
  degreeC = NULL, 
  degreeD = NULL,
  meanConversionToD = NULL, 
  meanConversionToC = NULL, 
  degreeLost = NULL,
  degreeLostC = NULL,
  degreeLostD = NULL
)


for(frac in c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)){
  #nodes in the top fractionCoop degrees will automatically be a cooperator
  fractionCoop = frac
  
    coopFrac = NULL
    avgCoop = NULL
    avgCoopFinal = NULL
    percentIsolation = NULL
    isolation = NULL
    percentIsolatedD = NULL
    nCommunities = NULL
    communitySize = NULL
    assortativityInitial = NULL
    assortativityFinal = NULL
    conversionRate = NULL
    conversionToD = NULL 
    conversionToC = NULL
    degreeC = NULL
    degreeD = NULL
    meanConversionToD = NULL 
    meanConversionToC = NULL
    degreeLost = NULL
    degreeLostC = NULL
    degreeLostD = NULL
    for(m in c(1:iterations)){
      # Section 1. NOTES, packages, and Parameters
      #Importing library
      library(igraph) # for network graphing
      library(reldist) # for gini calculatio
      library(boot) # for inv.logit calculation
      #Two prefixed functions
      #rank
      rank1 = function(x) {rank(x,na.last=NA,ties.method="average")[1]} #a smaller value has a smaller rank.
      #gini mean difference (a.k.a. mean difference: please refer to https://stat.ethz.ch/pipermail/r-help/2003-April/032782.html)
      gmd = function(x) {
       x1 = na.omit(x)
       n = length(x1)
      tmp = 0
       for (i in 1:n) {
       for (j in 1:n) {
       tmp <- tmp + abs(x1[i]-x1[j])
       }
       }
      answer = tmp/(n*n)
       return(answer)
      }
      # List of manipulating parameters of experiments
      #L : number of round
      #V : Visible or not
      #A : Income of a rich-group subject
      #B : Income of a poor-group subject
      #R : Probability to be assigned to a rich group
      #I : Number of the same-parameter trial
      #Example
      L = 10
      V = 0
      A = 700
      B = 300
      R = 0.5
      I = 0
      # List of fixed parameters of experiments (assumptions)
      #Rewiring rate = 0.3
      #GINI coefficient (can be known by A or B)
      GINI = 0*as.numeric(A==500) + 0.2*as.numeric(A %in% c(700,850)) + 0.4*as.numeric(A ==1150)
      #Collecting data frame (final output data frame)
      result = data.frame(round=0:L,n_par=NA,n_A=NA,avg_coop=NA,avg_degree=NA,avg_wealth=NA,gini=NA,gmd=NA,avg_coop_A=NA,avg_degree_A=NA,avg_wealth_A=NA,gini_A=NA,gmd_A=NA,avg_coop_B=NA,avg_degree_B=NA,avg_wealth_B=NA,gini_B=NA,gmd_B=NA,isolation=NA,percentIsolation=NA,meanConversionToD=NA,meanConversionToC=NA,degreeLost=NA,degreeLostC=NA,degreeLostD=NA)
      #_A is for a richer group and _B is for a poorer group
      
      
      #####################################################
      # Section 1.5: Practice rounds 1 to 2, to determine C/D in round 1
      
      N = 17 # median of the number of participants over rounds.
      node_rp0 = data.frame(ego_id=1:N, round=0)
      node_import = node_rp0
      
      for (k in 1:2){
        node_rX = node_import #Importing data
        node_rX$round = node_rX$round + 1
        node_rX[is.na(node_rX$prev_degree)==1,"prev_degree"] = 0
        node_rX[is.na(node_rX$prev_local_rate_coop)==1,"prev_local_rate_coop"] = 0
        #Only this calculation needs to change from Round 1
        if (k==1) {
          node_rX$prob_coop = inv.logit(1.099471) 
        } else {
          node_rX$prob_coop = inv.logit((-0.02339288) + (1.46068980)*as.numeric(node_rX$prev_coop==1))
        } 
        
        node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
        
        node_rX$prev_coop = node_rX$coop
          
        assign(paste("coop_rp",k, sep=""),node_rX$coop)
        
        #For the loop
        node_import = node_rX
      }
      
      #cooperation rate in the practice rounds
      coop_rp = apply(cbind(coop_rp1,coop_rp2),1,mean)

      #####################################################
      # Section 2: Round 0 (Agents and environments)
      #Node data generation
      N = 17 # median of the number of participants over rounds.
      node_r0 = data.frame(ego_id=1:N, round=0)
      node_r0$group = sample(c("rich","poor"),N,replace=TRUE,prob=c(R,1-R)) #R is defined as the probability to be assigned to the rich group
      node_r0$initial_wealth = ifelse(node_r0$group=="rich",A,B)
      #Link data generation
      ego_list = NULL
      for (i in 1:N) { ego_list = c(ego_list,rep(i,N)) }
      link_r0 = data.frame(ego_id=ego_list,alt_id=rep(1:N,N))
      link_r0 = link_r0[(link_r0$ego_id < link_r0$alt_id),] #The link was bidirectional, and thus the half and self are omitted.
      link_r0$connected = sample(0:1,dim(link_r0)[1],replace=TRUE,prob=c(0.7,0.3)) #Initial rewiring rate is fixed, 0.3
      link_r0c_ego = link_r0[link_r0$connected==1,]
      link_r0c_alt = link_r0[link_r0$connected==1,]
      colnames(link_r0c_alt) = c("alt_id","ego_id","connected")
      link_r0c = rbind(link_r0c_ego,link_r0c_alt) #this is bidirectional (double counted) for connected ties.
      link_r0c = link_r0c[order(link_r0c$ego_id),]
      link_r0c$alternumber = NA #putting the number for each alter in the same ego
      link_r0c[1,]$alternumber = 1
      for (i in 1:(dim(link_r0c)[1]-1))
        {if (link_r0c[i,]$ego_id == link_r0c[i+1,]$ego_id)
          {link_r0c[i+1,]$alternumber = link_r0c[i,]$alternumber + 1}
        else
          {link_r0c[i+1,]$alternumber = 1}
        #print(i)
        }
      link_r0c2 = reshape(link_r0c, direction = "wide", idvar=c("ego_id","connected"), timevar="alternumber")
      link_r0c2$initial_degree = apply(link_r0c2[,colnames(link_r0c2)[substr(colnames(link_r0c2),1,6) == "alt_id"]],1,function(x){length(na.omit(x))}) #Degree of each ego
      link_r0c2[is.na(link_r0c2$initial_degree)==1,"initial_degree"] = 0
      #Reflect the degree and initial local gini coefficient into the node data
      node_r0 = merge(x=node_r0,y=link_r0c2,all.x=TRUE,all.y=FALSE,by="ego_id")
      node_r0$initial_avg_env_wealth = NA
      node_r0$initial_local_gini = NA #local gini coefficient of the ego and connecting alters
      node_r0$initial_rel_rank = NA #local rank of ego among the ego and connecting alters (divided by the number of the go and connecting alters)
      for (i in 1:(dim(node_r0)[1])){
        node_r0[i,]$initial_avg_env_wealth = mean(na.omit(node_r0[node_r0$ego_id %in%
        node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6) %in% c("ego_id","alt_id")]],"initial_wealth"]))
        node_r0[i,]$initial_local_gini = gini(na.omit(node_r0[node_r0$ego_id %in% node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6)
        %in% c("ego_id","alt_id")]],"initial_wealth"]))
        node_r0[i,]$initial_rel_rank = rank1(na.omit(node_r0[node_r0$ego_id %in% node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6)
        %in% c("ego_id","alt_id")]],"initial_wealth"]))/length(na.omit(node_r0[node_r0$ego_id %in%
        node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6) %in% c("ego_id","alt_id")]],"initial_wealth"]))
        }
      #Finalization of round 0 and Visualization
      #plot(graph.data.frame(link_r0[link_r0$connected==1,],directed=F)) #plot.igraph
      node_r0$everIsolated = 0
      node_r0$maxDegreeLost = NA
      result[result$round==0,2:25] = c(length(node_r0$ego_id),length(node_r0[node_r0$group=="rich",]$ego_id),NA,mean(node_r0$initial_degree),mean(node_r0$initial_wealth),gini(node_r0$initial_wealth),gmd(node_r0$initial_wealth),NA,mean(node_r0[node_r0$group=="rich",]$initial_degree),mean(node_r0[node_r0$group=="rich",]$initial_wealth),gini(node_r0[node_r0$group=="rich",]$initial_wealth),gmd(node_r0[node_r0$group=="rich",]$initial_wealth),NA,mean(node_r0[node_r0$group=="poor",]$initial_degree),mean(node_r0[node_r0$group=="poor",]$initial_wealth),gini(node_r0[node_r0$group=="poor",]$initial_wealth),gmd(node_r0[node_r0$group=="poor",]$initial_wealth),
                                       as.numeric(ifelse(is.na(table(node_r0$initial_degree<=isolationDegree)["TRUE"]),0,1)),
                                       as.numeric(sum(node_r0$everIsolated)/length(node_r0$ego_id)),
                                       NA,
                                       NA,
                                       NA,NA,NA
                                       )
      
      #For the loop at the next round (for round 1, the initial one is the same as the previous [1 prior] one)
      node_import = node_r0
      node_import$initial_coop = NA
      node_import$prev_coop = NA
      node_import$prev_wealth = node_import$initial_wealth
      node_import$prev_degree = node_import$initial_degree
      node_import$prev_avg_env_wealth = node_import$initial_avg_env_wealth
      node_import$prev_local_gini = node_import$initial_local_gini
      node_import$prev_rel_rank = node_import$initial_rel_rank
      node_import$prev_local_rate_coop = NA
      link_import = link_r0
      
      
      #####################################################
      # Section 3: Rounds 1 to 10 or more (behaviors in simulation: the equation of cooperation is different at round 1 because of no history)
      #3-1: Cooperation phase
      for (k in 1:L)
      {
        node_rX = node_import #Importing data
        node_rX$round = node_rX$round + 1
        node_rX[is.na(node_rX$prev_degree)==1,"prev_degree"] = 0
        node_rX[is.na(node_rX$prev_local_rate_coop)==1,"prev_local_rate_coop"] = 0
        #Only this calculation needs to change from Round 1
        if (k==1) {
          node_rX$prob_coop = as.numeric(V==0)*inv.logit((-1.816665) + (2.086067)*coop_rp1 + (1.800153)*coop_rp2) + as.numeric(V==1)*inv.logit((-2.031577) + (2.427157)*coop_rp1 + (1.684193)*coop_rp2 + (-1.528851)*GINI)
          } else {
          node_rX$prob_coop = as.numeric(V==0 & node_rX$prev_coop==0)*inv.logit(-1.039916) + as.numeric(V==0 & node_rX$prev_coop==1)*inv.logit(2.062023) + as.numeric(V==1 & node_rX$prev_coop==0)*inv.logit((-0.2574838)*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-1.214198)*GINI + (2.508148)*GINI*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-0.9749075)) + as.numeric(V==1 & node_rX$prev_coop==1)*inv.logit((- 0.6197254)*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-0.7480261)*GINI + (1.169674)*GINI*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (1.356784))
          } 
        #####manipulate cooperation rate ar round 1 depending on the degree! (keep the total coopertion rate at round 1 constant)
        #####make it so that the nodes with the top fractionCoop*100 percentile of degrees will 100% be a cooperator, but the average percentage of being a cooperator will be the same
        if(k==1){
          prob_coop_df = NULL
          nodesCoop = NULL
          nodesCoop = node_rX$prev_degree>=quantile(node_rX$prev_degree,fractionCoop)
          prob_coop_df = 
            data.frame(
              prob_coop = rev(node_rX$prob_coop[order(coop_rp)]),
              node_number = c(which(nodesCoop),which(!nodesCoop))
            )
          node_rX$prob_coop = prob_coop_df[order(prob_coop_df$node_number),]$prob_coop
          node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
        } else {
          node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
        }
        if (k==1) {
          node_rX$initial_coop = node_rX$coop
          } else {
          node_rX$initial_coop = node_rX$initial_coop
          }
        node_rX$cost = (-50)*node_rX$coop*node_rX$prev_degree
        node_rX$n_coop_received = NA
        for (i in 1:(dim(node_rX)[1]))
          {
          node_rX[i,]$n_coop_received = sum(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) ==
          "alt_id"]],"coop"])
          }
        node_rX$benefit = 100*node_rX$n_coop_received
        node_rX$payoff = node_rX$cost + node_rX$benefit
        node_rX$wealth = node_rX$prev_wealth + node_rX$payoff
        node_rX$rel_rank = NA
        node_rX$local_rate_coop = NA
        for (i in 1:dim(node_rX)[1])
          {
          node_rX[i,]$rel_rank = rank1(na.omit(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in%
          c("ego_id","alt_id")]],"wealth"]))/length(na.omit(node_rX[node_rX$ego_id %in%
          node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
          node_rX[i,]$local_rate_coop = mean(na.omit(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in%
          c("ego_id","alt_id")]],"coop"]))
          }
        node_rX$growth = as.numeric((node_rX$wealth/node_rX$prev_wealth) > 1)
        node_rX = node_rX[,c("ego_id","round","group","prev_degree","initial_wealth","initial_local_gini","initial_coop","coop","wealth","rel_rank","local_rate_coop","growth","everIsolated","maxDegreeLost")] #Pruning the previous-round data (degree is not updating yet)
        
        #3-2: Rewiring phase
        # 30% of ties (unidirectional) are being rewired
        link_rX_1 = link_import #Importing data (bidirectioanl ego-alter [ego_id < alter_id])
        colnames(link_rX_1) = c("ego_id","alt_id","prev_connected")
        link_rX_1$challenge = sample(0:1,dim(link_rX_1)[1],replace=TRUE,prob=c(0.7,0.3)) # The bidirectional ties being rewired are selected (rewiring rate = 0.3).
        ego_node_data =
        node_rX[,c("ego_id","wealth","coop","prev_degree","initial_wealth","initial_local_gini","initial_coop","rel_rank","local_rate_coop","growth")]
        colnames(ego_node_data) =
        c("ego_id","ego_wealth","ego_coop","ego_prev_degree","ego_initial_wealth","ego_initial_local_gini","ego_initial_coop","ego_rel_rank","ego_local_rate_coop","ego_growth")
        alt_node_data =
        node_rX[,c("ego_id","wealth","coop","prev_degree","initial_wealth","initial_local_gini","initial_coop","rel_rank","local_rate_coop","growth")]
        colnames(alt_node_data) =
        c("alt_id","alt_wealth","alt_coop","alt_prev_degree","alt_initial_wealth","alt_initial_local_gini","alt_initial_coop","alt_rel_rank","alt_local_rate_coop","alt_growth")
        link_rX_2 = merge(x=link_rX_1,y=ego_node_data,all.x=TRUE,all.y=FALSE,by="ego_id")
        link_rX_3 = merge(x=link_rX_2,y=alt_node_data,all.x=TRUE,all.y=FALSE,by="alt_id")
        link_rX_3$choice = sample(c("ego","alt"),dim(link_rX_3)[1],replace=TRUE,prob=c(0.5,0.5)) #decision maker for breaking a link, which is a unilateral decision
        #ego_prob: probability of choosing to connect when challenged (asked)
        link_rX_3$ego_prob = inv.logit((0.5134401)*link_rX_3$prev_connected + (-0.852406)*link_rX_3$ego_coop + (2.96549)*link_rX_3$alt_coop + (-0.1808545))
        link_rX_3$alt_prob = inv.logit((0.5134401)*link_rX_3$prev_connected + (-0.852406)*link_rX_3$alt_coop + (2.96549)*link_rX_3$ego_coop + (-0.1808545))
        link_rX_3$prob_connect = ifelse(link_rX_3$prev_connected == 1, ifelse(link_rX_3$choice == "ego", link_rX_3$ego_prob,
        link_rX_3$alt_prob), link_rX_3$ego_prob*link_rX_3$alt_prob)
        link_rX_3$connect_update = apply(data.frame(link_rX_3$prob_connect),1, function(x) {sample(1:0,1,prob=c(x,(1-x)))})
        link_rX_3$connected = ifelse(link_rX_3$challenge==0,link_rX_3$prev_connected,link_rX_3$connect_update)
        link_rX = link_rX_3[,c("ego_id","alt_id","connected")] #pruning and data is updated
        #Reflect the degree and local gini coefficient into the node data
        link_rXc_ego = link_rX[link_rX$connected==1,]
        link_rXc_alt = link_rX[link_rX$connected==1,]
        colnames(link_rXc_alt) = c("alt_id","ego_id","connected")
        link_rXc = rbind(link_rXc_ego,link_rXc_alt)
        link_rXc = link_rXc[order(link_rXc$ego_id),]
        link_rXc$alternumber = NA
        link_rXc[1,]$alternumber = 1
        for (i in 1:(dim(link_rXc)[1]-1))
          {
            if (link_rXc[i,]$ego_id == link_rXc[i+1,]$ego_id)
            {
            link_rXc[i+1,]$alternumber = link_rXc[i,]$alternumber + 1
            }
            else
            {
            link_rXc[i+1,]$alternumber = 1
            }
            #print(i)
          }
        link_rXc2 = reshape(link_rXc, direction = "wide", idvar=c("ego_id","connected"), timevar="alternumber")
        link_rXc2$degree = apply(link_rXc2[,colnames(link_rXc2)[substr(colnames(link_rXc2),1,3) == "alt"]],1,function(x) {length(na.omit(x))})
        node_rX_final = merge(x=node_rX[,c("ego_id","round","group","initial_wealth","initial_local_gini","initial_coop","coop","wealth","growth","everIsolated","maxDegreeLost")],y=link_rXc2,all.x=TRUE,all.y=FALSE,by="ego_id")
        node_rX_final[is.na(node_rX_final$degree)==1,"degree"] = 0
        node_rX_final$avg_env_wealth = NA
        node_rX_final$local_gini = NA #needs to be updated because the social network changes at the rewiring phase
        node_rX_final$local_rate_coop = NA
        node_rX_final$rel_rank = NA
        for (i in 1:dim(node_rX_final)[1])
          {
          node_rX_final[i,]$avg_env_wealth = mean(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
          node_rX_final[i,]$local_gini = gini(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
          node_rX_final[i,]$local_rate_coop = mean(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"coop"]))
          node_rX_final[i,]$rel_rank = rank1(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in%
          c("ego_id","alt_id")]],"wealth"]))/length(na.omit(node_rX_final[node_rX_final$ego_id %in%
          node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
          node_rX_final[i,]$everIsolated = ifelse(node_rX_final[i,]$everIsolated==1,1,ifelse(node_rX_final[i,]$degree<=isolationDegree,1,0))
          node_rX_final[i,]$maxDegreeLost = pmax(node_r0[i,]$initial_degree - node_rX_final[i,]$degree, node_rX_final[i,]$maxDegreeLost, na.rm=TRUE)
        }
        
        
        #Finalization of round X and Visualization
        #plot(graph.data.frame(link_rX[link_rX$connected==1,],directed=F)) #plot.igraph
        result[result$round==k,2:25] =
        c(length(node_rX_final$ego_id),length(node_rX_final[node_rX_final$group=="rich",]$ego_id),mean(node_rX_final$coop),mean(node_rX_final$degree),mean(node_rX_final$wealth),gini(node_rX_final$wealth),gmd(node_rX_final$wealth),mean(node_rX_final[node_rX_final$group=="rich",]$coop),mean(node_rX_final[node_rX_final$group=="rich",]$degree),mean(node_rX_final[node_rX_final$group=="rich",]$wealth),gini(node_rX_final[node_rX_final$group=="rich",]$wealth),gmd(node_rX_final[node_rX_final$group=="rich",]$wealth),mean(node_rX_final[node_rX_final$group=="poor",]$coop),mean(node_rX_final[node_rX_final$group=="poor",]$degree),mean(node_rX_final[node_rX_final$group=="poor",]$wealth),gini(node_rX_final[node_rX_final$group=="poor",]$wealth),gmd(node_rX_final[node_rX_final$group=="poor",]$wealth),
                                       as.numeric(ifelse(is.na(table(node_rX_final$degree<=isolationDegree)["TRUE"]),0,1)),
                                       as.numeric(sum(node_rX_final$everIsolated)/length(node_rX_final$ego_id)),
          prop.table(table(node_rX_final[node_rX_final$initial_coop==1]$coop))["0"],
          prop.table(table(node_rX_final[node_rX_final$initial_coop==0]$coop))["1"],
          suppressWarnings({mean(node_rX_final$maxDegreeLost,na.rm=TRUE)}),
          suppressWarnings({mean(node_rX_final[node_rX_final$initial_coop==1]$maxDegreeLost,na.rm=TRUE)}),
          suppressWarnings({mean(node_rX_final[node_rX_final$initial_coop==0]$maxDegreeLost,na.rm=TRUE)})
          )
        
        #For the loop
        node_import = node_rX_final
        colnames(node_import)[colnames(node_import) %in%
        c("coop","wealth","growth","degree","avg_env_wealth","local_gini","local_rate_coop","rel_rank")] =
        c("prev_coop","prev_wealth","prev_growth","prev_degree","prev_avg_env_wealth","prev_local_gini","prev_local_rate_coop","prev_rel_rank")
        link_import = link_rX
        #print(paste0("Round ",k," is done."))
      }
      
      coopFrac[m] = fractionCoop
      avgCoop[m] = result[result$round==1,]$avg_coop
      avgCoopFinal[m] = result[result$round==10,]$avg_coop
      percentIsolation[m] = max(result[result$round>=1,]$percentIsolation)
      isolation[m] = max(result[result$round>=1,]$isolation)
      percentIsolatedD[m] = prop.table(table(node_rX_final[node_rX_final$everIsolated==1,]$initial_coop))["0"]
      
      
      link_rX_final = data.table::melt(setDT(node_rX_final),
                               measure = patterns('alt_id'),
                               variable.name = 'linkNumber', 
                               value.name = c('alt_id'))
      link_rX_final = data.frame(link_rX_final)[c("ego_id","alt_id")]
      link_rX_final = link_rX_final[complete.cases(link_rX_final),]
      link_rX_final = data.frame(t(unique(apply(link_rX_final, 1, function(x) sort(x))))) %>% distinct(X1, X2)
      node_g_final = data.frame(node_rX_final)[c("ego_id","initial_coop","coop")]
      node_g_final$initial_coop = factor(node_g_final$initial_coop)
      g_rX_final = graph_from_data_frame(link_rX_final, directed = FALSE, vertices=node_g_final)
      g_r0 = graph_from_data_frame(link_r0[link_r0$connected==1,][1:2], directed = FALSE, vertices=node_r0)
      nCommunities[m] = max(membership(cluster_louvain(g_rX_final)),na.rm=TRUE)
      communitySize[m] = mean(table(membership(cluster_louvain(g_rX_final))),na.rm=TRUE)
      assortativityInitial[m] = assortativity(g_r0, V(g_rX_final)$initial_coop == 1)
      assortativityFinal[m] = assortativity(g_rX_final, V(g_rX_final)$initial_coop == 1)
      conversionRate[m] = prop.table(table(V(g_rX_final)$coop == V(g_rX_final)$initial_coop))["FALSE"]
      conversionToD[m] = prop.table(table(V(g_rX_final)$coop[V(g_rX_final)$initial_coop==1]))["0"]
      conversionToC[m] = prop.table(table(V(g_rX_final)$coop[V(g_rX_final)$initial_coop==0]))["1"]
      degreeC[m] = mean(igraph::degree(g_rX_final)[V(g_rX_final)$coop==1],na.rm=TRUE)
      degreeD[m] = mean(igraph::degree(g_rX_final)[V(g_rX_final)$coop==0],na.rm=TRUE)
      meanConversionToD[m] = mean(result[result$round>=2,]$meanConversionToD, na.rm=TRUE)
      meanConversionToC[m] = mean(result[result$round>=2,]$meanConversionToC, na.rm=TRUE)
      degreeLost[m] = result[result$round==10,]$degreeLost
      degreeLostC[m] = result[result$round==10,]$degreeLostC
      degreeLostD[m] = result[result$round==10,]$degreeLostD
  
    }
    
df.netIntHighDegree = rbind(df.netIntHighDegree,
                  data.frame(
                    coopFrac = coopFrac,
                    avgCoop = avgCoop,
                    avgCoopFinal = avgCoopFinal,
                    percentIsolation = percentIsolation,
                    isolation = isolation,
                    percentIsolatedD = percentIsolatedD,
                    nCommunities = nCommunities,
                    communitySize = communitySize,
                    assortativityInitial = assortativityInitial,
                    assortativityFinal = assortativityFinal,
                    conversionRate = conversionRate,
                    conversionToD = conversionToD, 
                    conversionToC = conversionToC, 
                    degreeC = degreeC, 
                    degreeD = degreeD,
                    meanConversionToD = meanConversionToD, 
                    meanConversionToC = meanConversionToC,
                    degreeLost = degreeLost,
                    degreeLostC = degreeLostC,
                    degreeLostD = degreeLostD
                    ))
plot(g_r0,vertex.color=V(g_rX_final)$initial_coop,vertex.label=ifelse(is.na(V(g_rX_final)$initial_coop),"NA",ifelse(V(g_rX_final)$initial_coop==1,"C","D")),main=paste("fracCoop=",frac,", round 0",sep=""))
plot(g_rX_final,vertex.color=V(g_rX_final)$initial_coop,vertex.label=ifelse(is.na(V(g_rX_final)$initial_coop),"NA",ifelse(V(g_rX_final)$initial_coop==1,"C","D")),main=paste("fracCoop=",frac,", final round",sep=""))
}

sum.netIntHighDegree <- data.frame(
  df.netIntHighDegree %>% 
    group_by(coopFrac) %>%
      summarise(
        mean.isolation = mean(isolation),
        ci.isolation   = 1.96 * sd(isolation)/sqrt(n()),
        mean.percentIsolation = mean(percentIsolation),
        ci.percentIsolation   = 1.96 * sd(percentIsolation)/sqrt(n()),
        mean.percentIsolatedD = mean(percentIsolatedD,na.rm=TRUE),
        ci.percentIsolatedD   = 1.96 * sd(percentIsolatedD,na.rm=TRUE)/sqrt(sum(isolation)),
        mean.avgCoop = mean(avgCoop,na.rm=TRUE),
        ci.avgCoop   = 1.96 * sd(avgCoop,na.rm=TRUE)/sqrt(n()),
        mean.avgCoopFinal = mean(avgCoopFinal,na.rm=TRUE),
        ci.avgCoopFinal   = 1.96 * sd(avgCoopFinal,na.rm=TRUE)/sqrt(n()),
        mean.nCommunities = mean(nCommunities,na.rm=TRUE),
        ci.nCommunities   = 1.96 * sd(nCommunities,na.rm=TRUE)/sqrt(n()),
        mean.communitySize = mean(communitySize,na.rm=TRUE),
        ci.communitySize   = 1.96 * sd(communitySize,na.rm=TRUE)/sqrt(n()),
        mean.assortativityInitial = mean(assortativityInitial,na.rm=TRUE),
        ci.assortativityInitial   = 1.96 * sd(assortativityInitial,na.rm=TRUE)/sqrt(n()),
        mean.assortativityFinal = mean(assortativityFinal,na.rm=TRUE),
        ci.assortativityFinal   = 1.96 * sd(assortativityFinal,na.rm=TRUE)/sqrt(n()),
        mean.conversionRate = mean(conversionRate,na.rm=TRUE),
        ci.conversionRate   = 1.96 * sd(conversionRate,na.rm=TRUE)/sqrt(n()),
        mean.conversionToD = mean(conversionToD,na.rm=TRUE),
        ci.conversionToD   = 1.96 * sd(conversionToD,na.rm=TRUE)/sqrt(n()),
        mean.conversionToC = mean(conversionToC,na.rm=TRUE),
        ci.conversionToC   = 1.96 * sd(conversionToC,na.rm=TRUE)/sqrt(n()),
        mean.degreeC = mean(degreeC,na.rm=TRUE),
        ci.degreeC   = 1.96 * sd(degreeC,na.rm=TRUE)/sqrt(n()),
        mean.degreeD = mean(degreeD,na.rm=TRUE),
        ci.degreeD   = 1.96 * sd(degreeD,na.rm=TRUE)/sqrt(n()),
        mean.meanConversionToD = mean(meanConversionToD,na.rm=TRUE),
        ci.meanConversionToD   = 1.96 * sd(meanConversionToD,na.rm=TRUE)/sqrt(n()),
        mean.meanConversionToC = mean(meanConversionToC,na.rm=TRUE),
        ci.meanConversionToC   = 1.96 * sd(meanConversionToC,na.rm=TRUE)/sqrt(n()),
        mean.degreeLost = mean(degreeLost,na.rm=TRUE),
        ci.degreeLost   = 1.96 * sd(degreeLost,na.rm=TRUE)/sqrt(n()),
        mean.degreeLostC = mean(degreeLostC,na.rm=TRUE),
        ci.degreeLostC   = 1.96 * sd(degreeLostC,na.rm=TRUE)/sqrt(n()),
        mean.degreeLostD = mean(degreeLostD,na.rm=TRUE),
        ci.degreeLostD   = 1.96 * sd(degreeLostD,na.rm=TRUE)/sqrt(n())
        )
  )

kable(sum.netIntHighDegree[c(1:9)]) %>% kableExtra::kable_styling(font_size = 10)

coopFrac	mean.isolation	ci.isolation	mean.percentIsolation	ci.percentIsolation	mean.percentIsolatedD	ci.percentIsolatedD	mean.avgCoop	ci.avgCoop
0.0	0.554	0.0436142	0.0496471	0.0051501	0.8959716	0.0230048	0.7025882	0.0098537
0.1	0.626	0.0424550	0.0587059	0.0053488	0.9297577	0.0183035	0.7029412	0.0102805
0.2	0.710	0.0398138	0.0728235	0.0058606	0.9336816	0.0161372	0.6912941	0.0097417
0.3	0.686	0.0407223	0.0651765	0.0053147	0.9374613	0.0171653	0.7072941	0.0099539
0.4	0.680	0.0409294	0.0689412	0.0057807	0.9416667	0.0157113	0.7009412	0.0095685
0.5	0.612	0.0427560	0.0627059	0.0059694	0.9365223	0.0178315	0.7056471	0.0103070
0.6	0.632	0.0423144	0.0636471	0.0056710	0.9088578	0.0203857	0.7005882	0.0097025
0.7	0.570	0.0434388	0.0530588	0.0050686	0.9247967	0.0203289	0.7000000	0.0097519
0.8	0.612	0.0427560	0.0574118	0.0051598	0.9126031	0.0211305	0.6950588	0.0099950
0.9	0.590	0.0431543	0.0536471	0.0052337	0.9198639	0.0195573	0.6902353	0.0102957
1.0	0.576	0.0433611	0.0508235	0.0049005	0.8904454	0.0233704	0.6948235	0.0098804

kable(sum.netIntHighDegree[c(1,10:17)]) %>% kableExtra::kable_styling(font_size = 10) 

coopFrac	mean.avgCoopFinal	ci.avgCoopFinal	mean.nCommunities	ci.nCommunities	mean.communitySize	ci.communitySize	mean.assortativityInitial	ci.assortativityInitial
0.0	0.6977647	0.0096213	2.686	0.0463791	6.601667	0.1230765	-0.0745890	0.0108540
0.1	0.6940000	0.0100628	2.666	0.0446056	6.638500	0.1214974	-0.0744884	0.0112705
0.2	0.6880000	0.0101765	2.688	0.0463170	6.596000	0.1229046	-0.0613324	0.0105979
0.3	0.6988235	0.0101476	2.664	0.0463542	6.658333	0.1239272	-0.0698695	0.0105517
0.4	0.6911765	0.0099540	2.740	0.0447740	6.451500	0.1181238	-0.0631038	0.0113440
0.5	0.7040000	0.0099587	2.708	0.0463236	6.545000	0.1220174	-0.0628333	0.0121817
0.6	0.6944706	0.0096398	2.696	0.0460609	6.573333	0.1221946	-0.0487593	0.0117540
0.7	0.6962353	0.0096021	2.692	0.0455189	6.579000	0.1216156	-0.0736666	0.0113032
0.8	0.6940000	0.0102151	2.692	0.0484677	6.604500	0.1257684	-0.0599471	0.0120192
0.9	0.7031765	0.0099278	2.674	0.0457372	6.627167	0.1226940	-0.0685097	0.0112954
1.0	0.6976471	0.0097934	2.698	0.0449796	6.559167	0.1205845	-0.0629258	0.0111552

kable(sum.netIntHighDegree[c(1,18:25)]) %>% kableExtra::kable_styling(font_size = 10) 

coopFrac	mean.assortativityFinal	ci.assortativityFinal	mean.conversionRate	ci.conversionRate	mean.conversionToD	ci.conversionToD	mean.conversionToC	ci.conversionToC
0.0	-0.0578371	0.0053789	0.4130588	0.0100685	0.3042105	0.0112894	0.6994578	0.0180586
0.1	-0.0562541	0.0054400	0.4131765	0.0104958	0.3070120	0.0114594	0.6952051	0.0175166
0.2	-0.0575179	0.0060004	0.4214118	0.0101102	0.3141543	0.0111132	0.6812057	0.0185340
0.3	-0.0609540	0.0053862	0.4171765	0.0103065	0.3047862	0.0115100	0.7144351	0.0179432
0.4	-0.0570823	0.0056023	0.4140000	0.0105327	0.3101002	0.0116623	0.6855074	0.0175667
0.5	-0.0577792	0.0054489	0.4049412	0.0104693	0.2961342	0.0112179	0.7017002	0.0188656
0.6	-0.0512969	0.0052934	0.4178824	0.0107884	0.3082018	0.0113493	0.7027480	0.0181693
0.7	-0.0578332	0.0053707	0.4145882	0.0110138	0.3087364	0.0113233	0.6929049	0.0180900
0.8	-0.0568730	0.0056360	0.4116471	0.0100479	0.3024781	0.0110184	0.6889905	0.0187765
0.9	-0.0522759	0.0051748	0.4169412	0.0108513	0.2993473	0.0113404	0.7087871	0.0186078
1.0	-0.0552295	0.0052281	0.4091765	0.0109804	0.2996167	0.0114873	0.6919086	0.0184088

kable(sum.netIntHighDegree[c(1,26:33)]) %>% kableExtra::kable_styling(font_size = 10) 

coopFrac	mean.degreeC	ci.degreeC	mean.degreeD	ci.degreeD	mean.meanConversionToD	ci.meanConversionToD	mean.meanConversionToC	ci.meanConversionToC
0.0	11.08737	0.0860769	8.652930	0.1108792	0.2992567	0.0059296	0.6913703	0.0064344
0.1	11.16121	0.0819537	8.846247	0.1128566	0.2949080	0.0055140	0.6940359	0.0059650
0.2	11.07001	0.0874177	8.768482	0.1187724	0.2989252	0.0062004	0.6916606	0.0059999
0.3	11.17028	0.0812161	8.676856	0.1189345	0.2932942	0.0058708	0.6923120	0.0058214
0.4	11.06794	0.0822867	8.581765	0.1122660	0.3013115	0.0059873	0.6853333	0.0065409
0.5	11.07843	0.0836271	8.642444	0.1180098	0.2973419	0.0059277	0.6940483	0.0062217
0.6	11.02130	0.0808007	8.682728	0.1150781	0.3045915	0.0058418	0.6832126	0.0063607
0.7	11.06541	0.0855338	8.635409	0.1172260	0.3013643	0.0062912	0.6877214	0.0059058
0.8	11.10644	0.0836737	8.738735	0.1180338	0.3002972	0.0057637	0.6916027	0.0059537
0.9	11.10070	0.0798413	8.651301	0.1192510	0.2991218	0.0058637	0.6780647	0.0062944
1.0	11.11601	0.0830924	8.720908	0.1202779	0.2978674	0.0060136	0.6771242	0.0061063

kable(sum.netIntHighDegree[c(1,34:ncol(sum.netIntHighDegree))]) %>% kableExtra::kable_styling(font_size = 10) 

coopFrac	mean.degreeLost	ci.degreeLost	mean.degreeLostC	ci.degreeLostC	mean.degreeLostD	ci.degreeLostD
0.0	-0.8630441	0.0556564	-0.9342692	0.0535867	-0.6743237	0.0578146
0.1	-0.8798873	0.0528878	-0.9548275	0.0532751	-0.6890818	0.0481622
0.2	-0.7934505	0.0505385	-0.8179552	0.0524771	-0.7340351	0.0452436
0.3	-0.8715221	0.0507420	-0.9187971	0.0510982	-0.7410714	0.0481060
0.4	-0.8231250	0.0521859	-0.8737887	0.0513620	-0.7092962	0.0528086
0.5	-0.8631618	0.0510148	-0.8966292	0.0503789	-0.7698585	0.0520714
0.6	-0.8228456	0.0527547	-0.8881684	0.0513379	-0.6960424	0.0538364
0.7	-0.8568971	0.0546019	-0.9575688	0.0543539	-0.6666100	0.0512062
0.8	-0.8543471	0.0527343	-0.9387669	0.0541231	-0.7118319	0.0479064
0.9	-0.8261304	0.0522753	-0.9362643	0.0504786	-0.5984302	0.0503768
1.0	-0.8059966	0.0532102	-0.8965269	0.0515256	-0.6221925	0.0530861

compare_means(percentIsolation ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.              group1 group2             p       p.adj p.format p.signif method  
##    <chr>            <chr>  <chr>          <dbl>       <dbl> <chr>    <chr>    <chr>   
##  1 percentIsolation 0      0.1    0.00604       0.22        0.00604  **       Wilcoxon
##  2 percentIsolation 0      0.2    0.00000000134 0.000000073 1.3e-09  ****     Wilcoxon
##  3 percentIsolation 0      0.3    0.00000352    0.00018     3.5e-06  ****     Wilcoxon
##  4 percentIsolation 0      0.4    0.000000461   0.000024    4.6e-07  ****     Wilcoxon
##  5 percentIsolation 0      0.5    0.00329       0.13        0.00329  **       Wilcoxon
##  6 percentIsolation 0      0.6    0.000281      0.013       0.00028  ***      Wilcoxon
##  7 percentIsolation 0      0.7    0.243         1           0.24323  ns       Wilcoxon
##  8 percentIsolation 0      0.8    0.0135        0.46        0.01349  *        Wilcoxon
##  9 percentIsolation 0      0.9    0.218         1           0.21805  ns       Wilcoxon
## 10 percentIsolation 0      1      0.467         1           0.46666  ns       Wilcoxon
## # … with 45 more rows

compare_means(avgCoop ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.     group1 group2      p p.adj p.format p.signif method  
##    <chr>   <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 avgCoop 0      0.1    0.989      1 0.989    ns       Wilcoxon
##  2 avgCoop 0      0.2    0.0657     1 0.066    ns       Wilcoxon
##  3 avgCoop 0      0.3    0.713      1 0.713    ns       Wilcoxon
##  4 avgCoop 0      0.4    0.636      1 0.636    ns       Wilcoxon
##  5 avgCoop 0      0.5    0.740      1 0.740    ns       Wilcoxon
##  6 avgCoop 0      0.6    0.604      1 0.604    ns       Wilcoxon
##  7 avgCoop 0      0.7    0.489      1 0.489    ns       Wilcoxon
##  8 avgCoop 0      0.8    0.150      1 0.150    ns       Wilcoxon
##  9 avgCoop 0      0.9    0.0738     1 0.074    ns       Wilcoxon
## 10 avgCoop 0      1      0.160      1 0.160    ns       Wilcoxon
## # … with 45 more rows

compare_means(avgCoopFinal ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.          group1 group2     p p.adj p.format p.signif method  
##    <chr>        <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 avgCoopFinal 0      0.1    0.665     1 0.665    ns       Wilcoxon
##  2 avgCoopFinal 0      0.2    0.157     1 0.157    ns       Wilcoxon
##  3 avgCoopFinal 0      0.3    0.867     1 0.867    ns       Wilcoxon
##  4 avgCoopFinal 0      0.4    0.392     1 0.392    ns       Wilcoxon
##  5 avgCoopFinal 0      0.5    0.447     1 0.447    ns       Wilcoxon
##  6 avgCoopFinal 0      0.6    0.760     1 0.760    ns       Wilcoxon
##  7 avgCoopFinal 0      0.7    0.689     1 0.689    ns       Wilcoxon
##  8 avgCoopFinal 0      0.8    0.729     1 0.729    ns       Wilcoxon
##  9 avgCoopFinal 0      0.9    0.397     1 0.397    ns       Wilcoxon
## 10 avgCoopFinal 0      1      0.889     1 0.889    ns       Wilcoxon
## # … with 45 more rows

compare_means(nCommunities ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.          group1 group2      p p.adj p.format p.signif method  
##    <chr>        <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 nCommunities 0      0.1    0.628      1 0.628    ns       Wilcoxon
##  2 nCommunities 0      0.2    0.950      1 0.950    ns       Wilcoxon
##  3 nCommunities 0      0.3    0.519      1 0.519    ns       Wilcoxon
##  4 nCommunities 0      0.4    0.0903     1 0.090    ns       Wilcoxon
##  5 nCommunities 0      0.5    0.516      1 0.516    ns       Wilcoxon
##  6 nCommunities 0      0.6    0.754      1 0.754    ns       Wilcoxon
##  7 nCommunities 0      0.7    0.819      1 0.819    ns       Wilcoxon
##  8 nCommunities 0      0.8    0.964      1 0.964    ns       Wilcoxon
##  9 nCommunities 0      0.9    0.753      1 0.753    ns       Wilcoxon
## 10 nCommunities 0      1      0.660      1 0.660    ns       Wilcoxon
## # … with 45 more rows

compare_means(communitySize ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.           group1 group2      p p.adj p.format p.signif method  
##    <chr>         <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 communitySize 0      0.1    0.628      1 0.628    ns       Wilcoxon
##  2 communitySize 0      0.2    0.950      1 0.950    ns       Wilcoxon
##  3 communitySize 0      0.3    0.519      1 0.519    ns       Wilcoxon
##  4 communitySize 0      0.4    0.0903     1 0.090    ns       Wilcoxon
##  5 communitySize 0      0.5    0.516      1 0.516    ns       Wilcoxon
##  6 communitySize 0      0.6    0.754      1 0.754    ns       Wilcoxon
##  7 communitySize 0      0.7    0.819      1 0.819    ns       Wilcoxon
##  8 communitySize 0      0.8    0.964      1 0.964    ns       Wilcoxon
##  9 communitySize 0      0.9    0.753      1 0.753    ns       Wilcoxon
## 10 communitySize 0      1      0.660      1 0.660    ns       Wilcoxon
## # … with 45 more rows

compare_means(assortativityInitial ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.                  group1 group2       p p.adj p.format p.signif method  
##    <chr>                <chr>  <chr>    <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 assortativityInitial 0      0.1    0.921    1    0.9211   ns       Wilcoxon
##  2 assortativityInitial 0      0.2    0.129    1    0.1290   ns       Wilcoxon
##  3 assortativityInitial 0      0.3    0.739    1    0.7391   ns       Wilcoxon
##  4 assortativityInitial 0      0.4    0.309    1    0.3090   ns       Wilcoxon
##  5 assortativityInitial 0      0.5    0.447    1    0.4472   ns       Wilcoxon
##  6 assortativityInitial 0      0.6    0.00632  0.33 0.0063   **       Wilcoxon
##  7 assortativityInitial 0      0.7    0.914    1    0.9144   ns       Wilcoxon
##  8 assortativityInitial 0      0.8    0.178    1    0.1784   ns       Wilcoxon
##  9 assortativityInitial 0      0.9    0.472    1    0.4715   ns       Wilcoxon
## 10 assortativityInitial 0      1      0.234    1    0.2337   ns       Wilcoxon
## # … with 45 more rows

compare_means(assortativityFinal ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.                group1 group2      p p.adj p.format p.signif method  
##    <chr>              <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 assortativityFinal 0      0.1    0.481      1 0.4811   ns       Wilcoxon
##  2 assortativityFinal 0      0.2    0.984      1 0.9840   ns       Wilcoxon
##  3 assortativityFinal 0      0.3    0.596      1 0.5960   ns       Wilcoxon
##  4 assortativityFinal 0      0.4    0.533      1 0.5329   ns       Wilcoxon
##  5 assortativityFinal 0      0.5    0.754      1 0.7541   ns       Wilcoxon
##  6 assortativityFinal 0      0.6    0.0286     1 0.0286   *        Wilcoxon
##  7 assortativityFinal 0      0.7    0.785      1 0.7855   ns       Wilcoxon
##  8 assortativityFinal 0      0.8    0.416      1 0.4164   ns       Wilcoxon
##  9 assortativityFinal 0      0.9    0.0583     1 0.0583   ns       Wilcoxon
## 10 assortativityFinal 0      1      0.250      1 0.2504   ns       Wilcoxon
## # … with 45 more rows

compare_means(conversionRate ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.            group1 group2     p p.adj p.format p.signif method  
##    <chr>          <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 conversionRate 0      0.1    0.980     1 0.980    ns       Wilcoxon
##  2 conversionRate 0      0.2    0.237     1 0.237    ns       Wilcoxon
##  3 conversionRate 0      0.3    0.495     1 0.495    ns       Wilcoxon
##  4 conversionRate 0      0.4    0.944     1 0.944    ns       Wilcoxon
##  5 conversionRate 0      0.5    0.267     1 0.267    ns       Wilcoxon
##  6 conversionRate 0      0.6    0.546     1 0.546    ns       Wilcoxon
##  7 conversionRate 0      0.7    0.990     1 0.990    ns       Wilcoxon
##  8 conversionRate 0      0.8    0.951     1 0.951    ns       Wilcoxon
##  9 conversionRate 0      0.9    0.633     1 0.633    ns       Wilcoxon
## 10 conversionRate 0      1      0.546     1 0.546    ns       Wilcoxon
## # … with 45 more rows

compare_means(conversionToD ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.           group1 group2     p p.adj p.format p.signif method  
##    <chr>         <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 conversionToD 0      0.1    0.696     1 0.696    ns       Wilcoxon
##  2 conversionToD 0      0.2    0.180     1 0.180    ns       Wilcoxon
##  3 conversionToD 0      0.3    0.928     1 0.928    ns       Wilcoxon
##  4 conversionToD 0      0.4    0.546     1 0.546    ns       Wilcoxon
##  5 conversionToD 0      0.5    0.377     1 0.377    ns       Wilcoxon
##  6 conversionToD 0      0.6    0.588     1 0.588    ns       Wilcoxon
##  7 conversionToD 0      0.7    0.481     1 0.481    ns       Wilcoxon
##  8 conversionToD 0      0.8    0.843     1 0.843    ns       Wilcoxon
##  9 conversionToD 0      0.9    0.690     1 0.690    ns       Wilcoxon
## 10 conversionToD 0      1      0.462     1 0.462    ns       Wilcoxon
## # … with 45 more rows

compare_means(conversionToC ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.           group1 group2     p p.adj p.format p.signif method  
##    <chr>         <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 conversionToC 0      0.1    0.640     1 0.640    ns       Wilcoxon
##  2 conversionToC 0      0.2    0.245     1 0.245    ns       Wilcoxon
##  3 conversionToC 0      0.3    0.271     1 0.271    ns       Wilcoxon
##  4 conversionToC 0      0.4    0.290     1 0.290    ns       Wilcoxon
##  5 conversionToC 0      0.5    0.730     1 0.730    ns       Wilcoxon
##  6 conversionToC 0      0.6    0.779     1 0.779    ns       Wilcoxon
##  7 conversionToC 0      0.7    0.720     1 0.720    ns       Wilcoxon
##  8 conversionToC 0      0.8    0.518     1 0.518    ns       Wilcoxon
##  9 conversionToC 0      0.9    0.422     1 0.422    ns       Wilcoxon
## 10 conversionToC 0      1      0.588     1 0.588    ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeC ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.     group1 group2     p p.adj p.format p.signif method  
##    <chr>   <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeC 0      0.1    0.231     1 0.231    ns       Wilcoxon
##  2 degreeC 0      0.2    0.809     1 0.809    ns       Wilcoxon
##  3 degreeC 0      0.3    0.237     1 0.237    ns       Wilcoxon
##  4 degreeC 0      0.4    0.638     1 0.638    ns       Wilcoxon
##  5 degreeC 0      0.5    0.751     1 0.751    ns       Wilcoxon
##  6 degreeC 0      0.6    0.218     1 0.218    ns       Wilcoxon
##  7 degreeC 0      0.7    0.657     1 0.657    ns       Wilcoxon
##  8 degreeC 0      0.8    0.971     1 0.971    ns       Wilcoxon
##  9 degreeC 0      0.9    0.835     1 0.835    ns       Wilcoxon
## 10 degreeC 0      1      0.728     1 0.728    ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeD ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.     group1 group2      p p.adj p.format p.signif method  
##    <chr>   <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeD 0      0.1    0.0264     1 0.0264   *        Wilcoxon
##  2 degreeD 0      0.2    0.280      1 0.2802   ns       Wilcoxon
##  3 degreeD 0      0.3    0.674      1 0.6744   ns       Wilcoxon
##  4 degreeD 0      0.4    0.519      1 0.5189   ns       Wilcoxon
##  5 degreeD 0      0.5    0.694      1 0.6936   ns       Wilcoxon
##  6 degreeD 0      0.6    0.832      1 0.8318   ns       Wilcoxon
##  7 degreeD 0      0.7    0.898      1 0.8981   ns       Wilcoxon
##  8 degreeD 0      0.8    0.327      1 0.3266   ns       Wilcoxon
##  9 degreeD 0      0.9    0.847      1 0.8466   ns       Wilcoxon
## 10 degreeD 0      1      0.530      1 0.5304   ns       Wilcoxon
## # … with 45 more rows

compare_means(meanConversionToD ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.               group1 group2     p p.adj p.format p.signif method  
##    <chr>             <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 meanConversionToD 0      0.1    0.471     1 0.4714   ns       Wilcoxon
##  2 meanConversionToD 0      0.2    0.985     1 0.9851   ns       Wilcoxon
##  3 meanConversionToD 0      0.3    0.178     1 0.1780   ns       Wilcoxon
##  4 meanConversionToD 0      0.4    0.581     1 0.5805   ns       Wilcoxon
##  5 meanConversionToD 0      0.5    0.787     1 0.7868   ns       Wilcoxon
##  6 meanConversionToD 0      0.6    0.186     1 0.1861   ns       Wilcoxon
##  7 meanConversionToD 0      0.7    0.624     1 0.6237   ns       Wilcoxon
##  8 meanConversionToD 0      0.8    0.823     1 0.8227   ns       Wilcoxon
##  9 meanConversionToD 0      0.9    0.904     1 0.9043   ns       Wilcoxon
## 10 meanConversionToD 0      1      0.889     1 0.8886   ns       Wilcoxon
## # … with 45 more rows

compare_means(meanConversionToC ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.               group1 group2     p p.adj p.format p.signif method  
##    <chr>             <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 meanConversionToC 0      0.1    0.757     1 0.76     ns       Wilcoxon
##  2 meanConversionToC 0      0.2    0.780     1 0.78     ns       Wilcoxon
##  3 meanConversionToC 0      0.3    0.753     1 0.75     ns       Wilcoxon
##  4 meanConversionToC 0      0.4    0.709     1 0.71     ns       Wilcoxon
##  5 meanConversionToC 0      0.5    0.685     1 0.69     ns       Wilcoxon
##  6 meanConversionToC 0      0.6    0.693     1 0.69     ns       Wilcoxon
##  7 meanConversionToC 0      0.7    0.930     1 0.93     ns       Wilcoxon
##  8 meanConversionToC 0      0.8    0.717     1 0.72     ns       Wilcoxon
##  9 meanConversionToC 0      0.9    0.305     1 0.30     ns       Wilcoxon
## 10 meanConversionToC 0      1      0.409     1 0.41     ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeLost ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.        group1 group2      p p.adj p.format p.signif method  
##    <chr>      <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeLost 0      0.1    0.639      1 0.639    ns       Wilcoxon
##  2 degreeLost 0      0.2    0.0568     1 0.057    ns       Wilcoxon
##  3 degreeLost 0      0.3    0.860      1 0.860    ns       Wilcoxon
##  4 degreeLost 0      0.4    0.231      1 0.231    ns       Wilcoxon
##  5 degreeLost 0      0.5    0.773      1 0.773    ns       Wilcoxon
##  6 degreeLost 0      0.6    0.354      1 0.354    ns       Wilcoxon
##  7 degreeLost 0      0.7    0.965      1 0.965    ns       Wilcoxon
##  8 degreeLost 0      0.8    0.634      1 0.634    ns       Wilcoxon
##  9 degreeLost 0      0.9    0.388      1 0.388    ns       Wilcoxon
## 10 degreeLost 0      1      0.109      1 0.109    ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeLostC ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.         group1 group2       p p.adj p.format p.signif method  
##    <chr>       <chr>  <chr>    <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeLostC 0      0.1    0.518    1    0.5178   ns       Wilcoxon
##  2 degreeLostC 0      0.2    0.00942  0.49 0.0094   **       Wilcoxon
##  3 degreeLostC 0      0.3    0.726    1    0.7257   ns       Wilcoxon
##  4 degreeLostC 0      0.4    0.115    1    0.1155   ns       Wilcoxon
##  5 degreeLostC 0      0.5    0.200    1    0.2000   ns       Wilcoxon
##  6 degreeLostC 0      0.6    0.292    1    0.2918   ns       Wilcoxon
##  7 degreeLostC 0      0.7    0.541    1    0.5406   ns       Wilcoxon
##  8 degreeLostC 0      0.8    0.984    1    0.9843   ns       Wilcoxon
##  9 degreeLostC 0      0.9    0.860    1    0.8598   ns       Wilcoxon
## 10 degreeLostC 0      1      0.269    1    0.2691   ns       Wilcoxon
## # … with 45 more rows

compare_means(degreeLostD ~ coopFrac, data=df.netIntHighDegree)

## # A tibble: 55 × 8
##    .y.         group1 group2     p p.adj p.format p.signif method  
##    <chr>       <chr>  <chr>  <dbl> <dbl> <chr>    <chr>    <chr>   
##  1 degreeLostD 0      0.1    0.761     1 0.7607   ns       Wilcoxon
##  2 degreeLostD 0      0.2    0.338     1 0.3376   ns       Wilcoxon
##  3 degreeLostD 0      0.3    0.288     1 0.2881   ns       Wilcoxon
##  4 degreeLostD 0      0.4    0.566     1 0.5656   ns       Wilcoxon
##  5 degreeLostD 0      0.5    0.168     1 0.1685   ns       Wilcoxon
##  6 degreeLostD 0      0.6    0.573     1 0.5731   ns       Wilcoxon
##  7 degreeLostD 0      0.7    0.862     1 0.8622   ns       Wilcoxon
##  8 degreeLostD 0      0.8    0.495     1 0.4952   ns       Wilcoxon
##  9 degreeLostD 0      0.9    0.336     1 0.3358   ns       Wilcoxon
## 10 degreeLostD 0      1      0.592     1 0.5922   ns       Wilcoxon
## # … with 45 more rows

g.netIntHighDegree = ggbarplot(data=df.netIntHighDegree, x="coopFrac", y="percentIsolation", add = "mean_se") +
  stat_compare_means(ref.group = "0", label = "p.signif", label.y = 0.098, method="t.test") +  
  labs(
    title = "Isolation when cooperators are assigned to high-degree nodes",
    x = "Proportion of high-degree nodes assigned to cooperators",
    y = "Propoption of ever-isolated individuals") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(ylim=c(0,0.10))
g.netIntHighDegree